source: trunk/source/geometry/solids/CSG/src/G4Sphere.cc @ 836

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

import all except CVS
File size: 82.6 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: G4Sphere.cc,v 1.57.2.1 2008/04/23 09:05:23 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29//
30// class G4Sphere
31//
32// Implementation for G4Sphere class
33//
34// History:
35//
36// 22.07.05 O.Link    : Added check for intersection with double cone
37// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
38// 16.09.04 V.Grichine: bug fixed in SurfaceNormal(p), theta normals
39// 16.07.04 V.Grichine: bug fixed in DistanceToOut(p,v), Rmin go outside
40// 02.06.04 V.Grichine: bug fixed in DistanceToIn(p,v), on Rmax,Rmin go inside
41// 30.10.03 J.Apostolakis: new algorithm in Inside for SPhi-sections
42// 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi < SPhi, SPhi<0
43// 19.06.02 V.Grichine: bug fixed in Inside(p), && -> && fDTheta - kAngTolerance
44// 30.01.02 V.Grichine: bug fixed in Inside(p), && -> || at l.451
45// 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
46// 18.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
47// 25.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), phi intersections
48// 12.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), theta intersections
49// 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...)
50// 17.09.96 V.Grichine: final modifications to commit
51// 28.03.94 P.Kent: old C++ code converted to tolerant geometry
52// --------------------------------------------------------------------
53
54#include <assert.h>
55
56#include "G4Sphere.hh"
57
58#include "G4VoxelLimits.hh"
59#include "G4AffineTransform.hh"
60#include "G4GeometryTolerance.hh"
61
62#include "G4VPVParameterisation.hh"
63
64#include "Randomize.hh"
65
66#include "meshdefs.hh"
67
68#include "G4VGraphicsScene.hh"
69#include "G4VisExtent.hh"
70#include "G4Polyhedron.hh"
71#include "G4NURBS.hh"
72#include "G4NURBSbox.hh"
73
74using namespace CLHEP;
75
76// Private enum: Not for external use - used by distanceToOut
77
78enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTheta,kETheta};
79
80// used by normal
81
82enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSTheta,kNETheta};
83
84////////////////////////////////////////////////////////////////////////
85//
86// constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
87//             - note if pDPhi>2PI then reset to 2PI
88
89G4Sphere::G4Sphere( const G4String& pName,
90                          G4double pRmin, G4double pRmax,
91                          G4double pSPhi, G4double pDPhi,
92                          G4double pSTheta, G4double pDTheta )
93  : G4CSGSolid(pName)
94{
95  fEpsilon = 1.0e-14;
96
97  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
98  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
99
100  // Check radii
101
102  if (pRmin<pRmax&&pRmin>=0)
103  {
104    fRmin=pRmin; fRmax=pRmax;
105  }
106  else
107  {
108    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
109           << "        Invalide values for radii ! - "
110           << "        pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl;
111    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
112                "Invalid radii");
113  }
114
115  // Check angles
116
117  if (pDPhi>=twopi)
118  {
119    fDPhi=twopi;
120  }
121  else if (pDPhi>0)
122  {
123    fDPhi=pDPhi;
124  }
125  else
126  {
127    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
128           << "        Negative Z delta-Phi ! - "
129           << pDPhi << G4endl;
130    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
131                "Invalid DPhi.");
132  }
133
134  // Convert fSPhi to 0-2PI
135
136  if (pSPhi<0)
137  {
138    fSPhi=twopi-std::fmod(std::fabs(pSPhi),twopi);
139  }
140  else
141  {
142    fSPhi=std::fmod(pSPhi,twopi);
143  }
144
145  // Sphere is placed such that fSPhi+fDPhi>twopi !
146  // fSPhi could be < 0 !!?
147  //
148  if (fSPhi+fDPhi>twopi) fSPhi-=twopi;
149
150  // Check theta angles
151
152  if (pSTheta<0 || pSTheta>pi)
153  {
154    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl;
155    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
156                "stheta outside 0-PI range.");
157  }
158  else
159  {
160    fSTheta=pSTheta;
161  }
162
163  if (pDTheta+pSTheta>=pi)
164  {
165    fDTheta=pi-pSTheta;
166  }
167  else if (pDTheta>0)
168  {
169    fDTheta=pDTheta;
170  }
171  else
172  {
173    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
174           << "        Negative delta-Theta ! - "
175           << pDTheta << G4endl;
176    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
177                "Invalid pDTheta.");
178  }
179}
180
181///////////////////////////////////////////////////////////////////////
182//
183// Fake default constructor - sets only member data and allocates memory
184//                            for usage restricted to object persistency.
185//
186G4Sphere::G4Sphere( __void__& a )
187  : G4CSGSolid(a)
188{
189}
190
191/////////////////////////////////////////////////////////////////////
192//
193// Destructor
194
195G4Sphere::~G4Sphere()
196{
197}
198
199//////////////////////////////////////////////////////////////////////////
200//
201// Dispatch to parameterisation for replication mechanism dimension
202// computation & modification.
203
204void G4Sphere::ComputeDimensions(       G4VPVParameterisation* p,
205                                  const G4int n,
206                                  const G4VPhysicalVolume* pRep)
207{
208  p->ComputeDimensions(*this,n,pRep);
209}
210
211////////////////////////////////////////////////////////////////////////////
212//
213// Calculate extent under transform and specified limit
214
215G4bool G4Sphere::CalculateExtent( const EAxis pAxis,
216                                  const G4VoxelLimits& pVoxelLimit,
217                                  const G4AffineTransform& pTransform,
218                                        G4double& pMin, G4double& pMax ) const
219{
220  if ( fDPhi==twopi && fDTheta==pi)  // !pTransform.IsRotated() &&
221  {
222    // Special case handling for solid spheres-shells
223    // (rotation doesn't influence).
224    // Compute x/y/z mins and maxs for bounding box respecting limits,
225    // with early returns if outside limits. Then switch() on pAxis,
226    // and compute exact x and y limit for x/y case
227     
228    G4double xoffset,xMin,xMax;
229    G4double yoffset,yMin,yMax;
230    G4double zoffset,zMin,zMax;
231
232    G4double diff1,diff2,maxDiff,newMin,newMax;
233    G4double xoff1,xoff2,yoff1,yoff2;
234
235    xoffset=pTransform.NetTranslation().x();
236    xMin=xoffset-fRmax;
237    xMax=xoffset+fRmax;
238    if (pVoxelLimit.IsXLimited())
239    {
240      if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
241        || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
242      {
243        return false;
244      }
245      else
246      {
247        if (xMin<pVoxelLimit.GetMinXExtent())
248        {
249          xMin=pVoxelLimit.GetMinXExtent();
250        }
251        if (xMax>pVoxelLimit.GetMaxXExtent())
252        {
253          xMax=pVoxelLimit.GetMaxXExtent();
254        }
255      }
256    }
257
258    yoffset=pTransform.NetTranslation().y();
259    yMin=yoffset-fRmax;
260    yMax=yoffset+fRmax;
261    if (pVoxelLimit.IsYLimited())
262    {
263      if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
264        || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
265      {
266        return false;
267      }
268      else
269      {
270        if (yMin<pVoxelLimit.GetMinYExtent())
271        {
272          yMin=pVoxelLimit.GetMinYExtent();
273        }
274        if (yMax>pVoxelLimit.GetMaxYExtent())
275        {
276          yMax=pVoxelLimit.GetMaxYExtent();
277        }
278      }
279    }
280
281    zoffset=pTransform.NetTranslation().z();
282    zMin=zoffset-fRmax;
283    zMax=zoffset+fRmax;
284    if (pVoxelLimit.IsZLimited())
285    {
286      if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
287        || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
288      {
289        return false;
290      }
291      else
292      {
293        if (zMin<pVoxelLimit.GetMinZExtent())
294        {
295          zMin=pVoxelLimit.GetMinZExtent();
296        }
297        if (zMax>pVoxelLimit.GetMaxZExtent())
298        {
299          zMax=pVoxelLimit.GetMaxZExtent();
300        }
301      }
302    }
303
304    // Known to cut sphere
305
306    switch (pAxis)
307    {
308      case kXAxis:
309        yoff1=yoffset-yMin;
310        yoff2=yMax-yoffset;
311        if (yoff1>=0&&yoff2>=0)
312        {
313          // Y limits cross max/min x => no change
314          //
315          pMin=xMin;
316          pMax=xMax;
317        }
318        else
319        {
320          // Y limits don't cross max/min x => compute max delta x,
321          // hence new mins/maxs
322          //
323          diff1=std::sqrt(fRmax*fRmax-yoff1*yoff1);
324          diff2=std::sqrt(fRmax*fRmax-yoff2*yoff2);
325          maxDiff=(diff1>diff2) ? diff1:diff2;
326          newMin=xoffset-maxDiff;
327          newMax=xoffset+maxDiff;
328          pMin=(newMin<xMin) ? xMin : newMin;
329          pMax=(newMax>xMax) ? xMax : newMax;
330        }
331        break;
332      case kYAxis:
333        xoff1=xoffset-xMin;
334        xoff2=xMax-xoffset;
335        if (xoff1>=0&&xoff2>=0)
336        {
337          // X limits cross max/min y => no change
338          //
339          pMin=yMin;
340          pMax=yMax;
341        }
342        else
343        {
344          // X limits don't cross max/min y => compute max delta y,
345          // hence new mins/maxs
346          //
347          diff1=std::sqrt(fRmax*fRmax-xoff1*xoff1);
348          diff2=std::sqrt(fRmax*fRmax-xoff2*xoff2);
349          maxDiff=(diff1>diff2) ? diff1:diff2;
350          newMin=yoffset-maxDiff;
351          newMax=yoffset+maxDiff;
352          pMin=(newMin<yMin) ? yMin : newMin;
353          pMax=(newMax>yMax) ? yMax : newMax;
354        }
355        break;
356      case kZAxis:
357        pMin=zMin;
358        pMax=zMax;
359        break;
360      default:
361        break;
362    }
363    pMin-=kCarTolerance;
364    pMax+=kCarTolerance;
365
366    return true; 
367  }
368  else       // Transformed cutted sphere
369  {
370    G4int i,j,noEntries,noBetweenSections;
371    G4bool existsAfterClip=false;
372
373    // Calculate rotated vertex coordinates
374
375    G4ThreeVectorList* vertices;
376    G4int  noPolygonVertices ;
377    vertices=CreateRotatedVertices(pTransform,noPolygonVertices);
378
379    pMin=+kInfinity;
380    pMax=-kInfinity;
381
382    noEntries=vertices->size();  // noPolygonVertices*noPhiCrossSections
383    noBetweenSections=noEntries-noPolygonVertices;
384
385    G4ThreeVectorList ThetaPolygon ;
386    for (i=0;i<noEntries;i+=noPolygonVertices)
387    {
388      for(j=0;j<(noPolygonVertices/2)-1;j++)
389      {
390        ThetaPolygon.push_back((*vertices)[i+j]) ;     
391        ThetaPolygon.push_back((*vertices)[i+j+1]) ;     
392        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;     
393        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;     
394        CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
395        ThetaPolygon.clear() ;
396      }
397    }
398    for (i=0;i<noBetweenSections;i+=noPolygonVertices)
399    {
400      for(j=0;j<noPolygonVertices-1;j++)
401      {
402        ThetaPolygon.push_back((*vertices)[i+j]) ;     
403        ThetaPolygon.push_back((*vertices)[i+j+1]) ;     
404        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;     
405        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;     
406        CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
407        ThetaPolygon.clear() ;
408      }
409      ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;     
410      ThetaPolygon.push_back((*vertices)[i]) ; 
411      ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;     
412      ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;     
413      CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
414      ThetaPolygon.clear() ;
415    }
416     
417    if (pMin!=kInfinity || pMax!=-kInfinity)
418    {
419      existsAfterClip=true;
420
421      // Add 2*tolerance to avoid precision troubles
422      //
423      pMin-=kCarTolerance;
424      pMax+=kCarTolerance;
425    }
426    else
427    {
428      // Check for case where completely enveloping clipping volume
429      // If point inside then we are confident that the solid completely
430      // envelopes the clipping volume. Hence set min/max extents according
431      // to clipping volume extents along the specified axis.
432
433      G4ThreeVector clipCentre(
434          (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
435          (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
436          (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
437       
438      if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
439      {
440        existsAfterClip=true;
441        pMin=pVoxelLimit.GetMinExtent(pAxis);
442        pMax=pVoxelLimit.GetMaxExtent(pAxis);
443      }
444    }
445    delete vertices;
446    return existsAfterClip;
447  }
448}
449
450///////////////////////////////////////////////////////////////////////////
451//
452// Return whether point inside/outside/on surface
453// Split into radius, phi, theta checks
454// Each check modifies `in', or returns as approprate
455
456EInside G4Sphere::Inside( const G4ThreeVector& p ) const
457{
458  G4double rho,rho2,rad2,tolRMin,tolRMax;
459  G4double pPhi,pTheta;
460  EInside in=kOutside;
461
462  rho2 = p.x()*p.x() + p.y()*p.y() ;
463  rad2 = rho2 + p.z()*p.z() ;
464
465  //  if(rad2 >= 1.369e+19) DBG();
466  //  G4double rad = std::sqrt(rad2);
467  // Check radial surfaces
468  // sets `in'
469
470  if ( fRmin ) tolRMin = fRmin + kRadTolerance*0.5;
471  else         tolRMin = 0 ;
472 
473  tolRMax = fRmax - kRadTolerance*0.5 ;
474  //  const G4double  fractionTolerance = 1.0e-12;
475  const G4double  flexRadMaxTolerance = // kRadTolerance;
476    std::max(kRadTolerance, fEpsilon * fRmax);
477
478  const G4double  Rmax_minus = fRmax - flexRadMaxTolerance*0.5;
479  const G4double  flexRadMinTolerance = std::max(kRadTolerance, 
480                     fEpsilon * fRmin);
481  const G4double  Rmin_plus = (fRmin > 0) ? fRmin + flexRadMinTolerance*0.5 : 0 ;
482   
483if(rad2 <= Rmax_minus*Rmax_minus && rad2 >= Rmin_plus*Rmin_plus) in = kInside ; 
484
485// if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kInside ;
486  // if ( rad <= tolRMax && rad >= tolRMin )  in = kInside ;
487  else
488  {
489    tolRMax = fRmax + kRadTolerance*0.5 ;
490    tolRMin = fRmin - kRadTolerance*0.5 ;
491
492    if ( tolRMin < 0.0 ) tolRMin = 0.0 ;
493   
494     if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kSurface ;
495    //  if ( rad <= tolRMax && rad >= tolRMin )  in = kSurface ;
496    else                                                return in = kOutside ;
497  }
498
499  // Phi boundaries   : Do not check if it has no phi boundary!
500  // (in != kOutside). It is new J.Apostolakis proposal of 30.10.03
501
502  if ( ( fDPhi < twopi - kAngTolerance ) &&
503       ( (p.x() != 0.0 ) || (p.y() != 0.0) ) )
504  {
505    pPhi = std::atan2(p.y(),p.x()) ;
506
507    if      ( pPhi < fSPhi - kAngTolerance*0.5  )         pPhi += twopi ; 
508    else if ( pPhi > fSPhi + fDPhi + kAngTolerance*0.5 )  pPhi -= twopi; 
509   
510    if ((pPhi < fSPhi - kAngTolerance*0.5) || 
511        (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) )  return in = kOutside ;
512   
513    else if (in == kInside)  // else it's kSurface anyway already
514    {
515      if ( (pPhi < fSPhi + kAngTolerance*0.5) || 
516           (pPhi > fSPhi + fDPhi - kAngTolerance*0.5) )      in = kSurface ;       
517    }
518  }
519
520  // Theta bondaries
521  // (in!=kOutside)
522 
523  if ( (rho2 || p.z()) && fDTheta < pi - kAngTolerance*0.5 )
524  {
525    rho    = std::sqrt(rho2);
526    pTheta = std::atan2(rho,p.z());
527
528    if ( in == kInside )
529    {
530      if ( (pTheta < fSTheta + kAngTolerance*0.5)
531        || (pTheta > fSTheta + fDTheta - kAngTolerance*0.5) )
532      {
533        if ( (pTheta >= fSTheta - kAngTolerance*0.5)
534          && (pTheta <= fSTheta + fDTheta + kAngTolerance*0.5) )
535        {
536          in = kSurface ;
537        }
538        else
539        {
540          in = kOutside ;
541        }
542      }
543    }
544    else
545    {
546      if ( (pTheta < fSTheta - kAngTolerance*0.5)
547        || (pTheta > fSTheta + fDTheta + kAngTolerance*0.5) )
548      {
549        in = kOutside ;
550      }
551    }
552  }
553  return in;
554}
555
556/////////////////////////////////////////////////////////////////////
557//
558// Return unit normal of surface closest to p
559// - note if point on z axis, ignore phi divided sides
560// - unsafe if point close to z axis a rmin=0 - no explicit checks
561
562G4ThreeVector G4Sphere::SurfaceNormal( const G4ThreeVector& p ) const
563{
564  G4int noSurfaces = 0; 
565  G4double rho, rho2, rad, pTheta, pPhi=0.;
566  G4double distRMin = kInfinity;
567  G4double distSPhi = kInfinity, distEPhi = kInfinity;
568  G4double distSTheta = kInfinity, distETheta = kInfinity;
569  G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
570  G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
571  G4ThreeVector norm, sumnorm(0.,0.,0.);
572
573  rho2 = p.x()*p.x()+p.y()*p.y();
574  rad  = std::sqrt(rho2+p.z()*p.z());
575  rho  = std::sqrt(rho2);
576
577  G4double    distRMax = std::fabs(rad-fRmax);
578  if (fRmin)  distRMin = std::fabs(rad-fRmin);
579   
580  if ( rho && (fDPhi < twopi || fDTheta < pi) )
581  {
582    pPhi = std::atan2(p.y(),p.x());
583
584    if(pPhi  < fSPhi-dAngle)           pPhi     += twopi;
585    else if(pPhi > fSPhi+fDPhi+dAngle) pPhi     -= twopi;
586  }
587  if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
588  {
589    if ( rho )
590    {
591      distSPhi = std::fabs( pPhi - fSPhi ); 
592      distEPhi = std::fabs(pPhi-fSPhi-fDPhi); 
593    }
594    else if( !fRmin )
595    {
596      distSPhi = 0.; 
597      distEPhi = 0.; 
598    }
599    nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
600    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
601  }       
602  if ( fDTheta < pi ) // && rad ) // old limitation against (0,0,0)
603  {
604    if ( rho )
605    {
606      pTheta     = std::atan2(rho,p.z());
607      distSTheta = std::fabs(pTheta-fSTheta); 
608      distETheta = std::fabs(pTheta-fSTheta-fDTheta);
609 
610      nTs = G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi),
611                        -std::cos(fSTheta)*std::sin(pPhi),
612                         std::sin(fSTheta)               );
613      nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi),
614                         std::cos(fSTheta+fDTheta)*std::sin(pPhi),
615                        -std::sin(fSTheta+fDTheta)               );   
616    }
617    else if( !fRmin )
618    {
619      if ( fSTheta )                distSTheta = 0.;
620      if ( fSTheta + fDTheta < pi ) distETheta = 0.;
621    }   
622  }
623  if( rad )  nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad);
624
625  if( distRMax <= delta )
626  {
627    noSurfaces ++;
628    sumnorm += nR;
629  }
630  if( fRmin && distRMin <= delta )
631  {
632    noSurfaces ++;
633    sumnorm -= nR;
634  }
635  if( fDPhi < twopi )   
636  {
637    if (distSPhi <= dAngle)
638    {
639      noSurfaces ++;
640      sumnorm += nPs;
641    }
642    if (distEPhi <= dAngle) 
643    {
644      noSurfaces ++;
645      sumnorm += nPe;
646    }
647  }
648  if ( fDTheta < pi )
649  {
650    if (distSTheta <= dAngle && fSTheta > 0.)
651    {
652      noSurfaces ++;
653      if( rad <= delta && fDPhi >= twopi) sumnorm += nZ;
654      else                                sumnorm += nTs;
655    }
656    if (distETheta <= dAngle && fSTheta+fDTheta < pi) 
657    {
658      noSurfaces ++;
659      if( rad <= delta && fDPhi >= twopi) sumnorm -= nZ;
660      else                                sumnorm += nTe;
661      if(sumnorm.z() == 0.)               sumnorm += nZ;
662    }
663  }
664  if ( noSurfaces == 0 )
665  {
666#ifdef G4CSGDEBUG
667    G4Exception("G4Sphere::SurfaceNormal(p)", "Notification", JustWarning, 
668                "Point p is not on surface !?" ); 
669#endif
670     norm = ApproxSurfaceNormal(p);
671  }
672  else if ( noSurfaces == 1 ) norm = sumnorm;
673  else                        norm = sumnorm.unit();
674  return norm;
675}
676
677
678/////////////////////////////////////////////////////////////////////////////////////////////
679//
680// Algorithm for SurfaceNormal() following the original specification
681// for points not on the surface
682
683G4ThreeVector G4Sphere::ApproxSurfaceNormal( const G4ThreeVector& p ) const
684{
685  ENorm side;
686  G4ThreeVector norm;
687  G4double rho,rho2,rad,pPhi,pTheta;
688  G4double distRMin,distRMax,distSPhi,distEPhi,
689           distSTheta,distETheta,distMin;
690
691  rho2=p.x()*p.x()+p.y()*p.y();
692  rad=std::sqrt(rho2+p.z()*p.z());
693  rho=std::sqrt(rho2);
694
695  //
696  // Distance to r shells
697  //
698
699  distRMax=std::fabs(rad-fRmax);
700  if (fRmin)
701  {
702    distRMin=std::fabs(rad-fRmin);
703     
704    if (distRMin<distRMax)
705    {
706      distMin=distRMin;
707      side=kNRMin;
708    }
709    else
710    {   
711      distMin=distRMax;
712      side=kNRMax;
713    }
714  }
715  else
716  {
717    distMin=distRMax;
718    side=kNRMax;
719  }
720
721  //
722  // Distance to phi planes
723  //
724  // Protected against (0,0,z)
725   
726  pPhi = std::atan2(p.y(),p.x());
727  if (pPhi<0) pPhi += twopi;
728
729  if (fDPhi<twopi&&rho)
730  {
731    if (fSPhi<0)
732    {
733      distSPhi=std::fabs(pPhi-(fSPhi+twopi))*rho;
734    }
735    else
736    {
737      distSPhi=std::fabs(pPhi-fSPhi)*rho;
738    }
739
740    distEPhi=std::fabs(pPhi-fSPhi-fDPhi)*rho;
741
742    // Find new minimum
743    //
744    if (distSPhi<distEPhi)
745    {
746      if (distSPhi<distMin)
747      {
748        distMin=distSPhi;
749        side=kNSPhi;
750      }
751    }
752    else
753    {
754      if (distEPhi<distMin)
755      {
756        distMin=distEPhi;
757        side=kNEPhi;
758      }
759    }
760  }
761
762  //
763  // Distance to theta planes
764  //
765
766  if (fDTheta<pi&&rad)
767  {
768    pTheta=std::atan2(rho,p.z());
769    distSTheta=std::fabs(pTheta-fSTheta)*rad;
770    distETheta=std::fabs(pTheta-fSTheta-fDTheta)*rad;
771
772    // Find new minimum
773    //
774    if (distSTheta<distETheta)
775    {
776      if (distSTheta<distMin)
777      {
778        distMin = distSTheta ;
779        side = kNSTheta ;
780      }
781    }
782    else
783    {
784      if (distETheta<distMin)
785      {
786        distMin = distETheta ;
787        side = kNETheta ;
788      }
789    }
790  }
791
792  switch (side)
793  {
794    case kNRMin:      // Inner radius
795      norm=G4ThreeVector(-p.x()/rad,-p.y()/rad,-p.z()/rad);
796      break;
797    case kNRMax:      // Outer radius
798      norm=G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad);
799      break;
800    case kNSPhi:
801      norm=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
802      break;
803    case kNEPhi:
804      norm=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
805      break;
806    case kNSTheta:
807      norm=G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi),
808                         -std::cos(fSTheta)*std::sin(pPhi),
809                          std::sin(fSTheta)            );
810      //  G4cout<<G4endl<<" case kNSTheta:"<<G4endl;
811      //  G4cout<<"pPhi = "<<pPhi<<G4endl;
812      //  G4cout<<"rad  = "<<rad<<G4endl;
813      //  G4cout<<"pho  = "<<rho<<G4endl;
814      //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
815      //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
816      break;
817    case kNETheta:
818      norm=G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi),
819                          std::cos(fSTheta+fDTheta)*std::sin(pPhi),
820                         -std::sin(fSTheta+fDTheta)              );
821
822      //  G4cout<<G4endl<<" case kNETheta:"<<G4endl;
823      //  G4cout<<"pPhi = "<<pPhi<<G4endl;
824      //  G4cout<<"rad  = "<<rad<<G4endl;
825      //  G4cout<<"pho  = "<<rho<<G4endl;
826      //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
827      //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
828      break;
829    default:
830      DumpInfo();
831      G4Exception("G4Sphere::ApproxSurfaceNormal()", "Notification", JustWarning,
832                  "Undefined side for valid surface normal to solid.");
833      break;   
834  } // end case
835
836  return norm;
837}
838
839///////////////////////////////////////////////////////////////////////////////
840//
841// Calculate distance to shape from outside, along normalised vector
842// - return kInfinity if no intersection, or intersection distance <= tolerance
843//
844// -> If point is outside outer radius, compute intersection with rmax
845//        - if no intersection return
846//        - if  valid phi,theta return intersection Dist
847//
848// -> If shell, compute intersection with inner radius, taking largest +ve root
849//        - if valid phi,theta, save intersection
850//
851// -> If phi segmented, compute intersection with phi half planes
852//        - if valid intersection(r,theta), return smallest intersection of
853//          inner shell & phi intersection
854//
855// -> If theta segmented, compute intersection with theta cones
856//        - if valid intersection(r,phi), return smallest intersection of
857//          inner shell & theta intersection
858//
859//
860// NOTE:
861// - `if valid' (above) implies tolerant checking of intersection points
862//
863// OPT:
864// Move tolIO/ORmin/RMax2 precalcs to where they are needed -
865// not required for most cases.
866// Avoid atan2 for non theta cut G4Sphere.
867
868G4double G4Sphere::DistanceToIn( const G4ThreeVector& p,
869                                 const G4ThreeVector& v  ) const
870{
871  G4double snxt = kInfinity ;      // snxt = default return value
872
873  G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
874
875  G4double tolIRMin2, tolORMin2, tolORMax2, tolIRMax2 ;
876  G4double tolSTheta=0., tolETheta=0. ;
877
878  // Intersection point
879
880  G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
881
882  // Phi intersection
883
884  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi , Comp ; 
885
886  // Phi flag and precalcs
887
888  G4bool segPhi ;       
889  G4double hDPhi, hDPhiOT, hDPhiIT, cPhi, sinCPhi=0., cosCPhi=0. ; 
890  G4double cosHDPhiOT=0., cosHDPhiIT=0. ;
891  G4double Dist, cosPsi ;
892
893  G4bool segTheta ;                             // Theta flag and precals
894  G4double tanSTheta, tanETheta ;
895  G4double tanSTheta2, tanETheta2 ;
896  G4double dist2STheta, dist2ETheta ;
897  G4double t1, t2, b, c, d2, d, s = kInfinity ;
898
899  // General Precalcs
900
901  rho2 = p.x()*p.x() + p.y()*p.y() ;
902  rad2 = rho2 + p.z()*p.z() ;
903  pTheta = std::atan2(std::sqrt(rho2),p.z()) ;
904
905  pDotV2d = p.x()*v.x() + p.y()*v.y() ;
906  pDotV3d = pDotV2d + p.z()*v.z() ;
907
908  // Radial Precalcs
909
910  if (fRmin > kRadTolerance*0.5)
911  {
912    tolORMin2=(fRmin-kRadTolerance*0.5)*(fRmin-kRadTolerance*0.5);
913  }
914  else
915  {
916    tolORMin2 = 0 ;
917  }
918  tolIRMin2 = (fRmin+kRadTolerance*0.5)*(fRmin+kRadTolerance*0.5) ;
919  tolORMax2 = (fRmax+kRadTolerance*0.5)*(fRmax+kRadTolerance*0.5) ;
920  tolIRMax2 = (fRmax-kRadTolerance*0.5)*(fRmax-kRadTolerance*0.5) ;
921
922  // Set phi divided flag and precalcs
923
924  if (fDPhi < twopi)
925  {
926    segPhi = true ;
927    hDPhi = 0.5*fDPhi ;    // half delta phi
928    cPhi = fSPhi + hDPhi ;
929
930    hDPhiOT = hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi
931    hDPhiIT = hDPhi-0.5*kAngTolerance;
932
933    sinCPhi    = std::sin(cPhi) ;
934    cosCPhi    = std::cos(cPhi) ;
935    cosHDPhiOT = std::cos(hDPhiOT) ;
936    cosHDPhiIT = std::cos(hDPhiIT) ;
937  }
938  else
939  {
940    segPhi = false ;
941  }
942
943  // Theta precalcs
944   
945  if (fDTheta < pi )
946  {
947    segTheta  = true ;
948    tolSTheta = fSTheta - kAngTolerance*0.5 ;
949    tolETheta = fSTheta + fDTheta + kAngTolerance*0.5 ;
950  }
951  else
952  {
953    segTheta = false ;
954  }
955
956  // Outer spherical shell intersection
957  // - Only if outside tolerant fRmax
958  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
959  // - No intersect -> no intersection with G4Sphere
960  //
961  // Shell eqn: x^2+y^2+z^2=RSPH^2
962  //
963  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
964  //
965  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
966  // =>      rad2        +2s(pDotV3d)       +s^2                =R^2
967  //
968  // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
969
970  c = rad2 - fRmax*fRmax ;
971  const G4double  flexRadMaxTolerance = // kRadTolerance;
972    std::max(kRadTolerance, fEpsilon * fRmax);
973
974  //  if (c > kRadTolerance*fRmax)
975  if (c > flexRadMaxTolerance*fRmax)
976  {
977    // If outside toleranct boundary of outer G4Sphere
978    // [should be std::sqrt(rad2)-fRmax > kRadTolerance*0.5]
979
980    d2 = pDotV3d*pDotV3d - c ;
981
982    if ( d2 >= 0 )
983    {
984      s = -pDotV3d - std::sqrt(d2) ;
985
986      if (s >= 0 )
987      {
988        xi   = p.x() + s*v.x() ;
989        yi   = p.y() + s*v.y() ;
990        rhoi = std::sqrt(xi*xi + yi*yi) ;
991
992        if (segPhi && rhoi)    // Check phi intersection
993        {
994          cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
995
996          if (cosPsi >= cosHDPhiOT)
997          {
998            if (segTheta)   // Check theta intersection
999            {
1000              zi = p.z() + s*v.z() ;
1001
1002              // rhoi & zi can never both be 0
1003              // (=>intersect at origin =>fRmax=0)
1004              //
1005              iTheta = std::atan2(rhoi,zi) ;
1006              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1007              {
1008                return snxt = s ;
1009              }
1010            }
1011            else
1012            {
1013              return snxt=s;
1014            }
1015          }
1016        }
1017        else
1018        {
1019          if (segTheta)    // Check theta intersection
1020          {
1021            zi = p.z() + s*v.z() ;
1022
1023            // rhoi & zi can never both be 0
1024            // (=>intersect at origin => fRmax=0 !)
1025            //
1026            iTheta = std::atan2(rhoi,zi) ;
1027            if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1028            {
1029              return snxt=s;
1030            }
1031          }
1032          else
1033          {
1034            return snxt = s ;
1035          }
1036        }         
1037      }
1038    }
1039    else    // No intersection with G4Sphere
1040    {
1041      return snxt=kInfinity;
1042    }
1043  }
1044  else
1045  {
1046    // Inside outer radius
1047    // check not inside, and heading through G4Sphere (-> 0 to in)
1048
1049    d2 = pDotV3d*pDotV3d - c ;
1050
1051    // if (rad2 > tolIRMin2 && pDotV3d < 0 )
1052
1053    if (rad2 > tolIRMax2 && ( d2 >= flexRadMaxTolerance*fRmax && pDotV3d < 0 ) )
1054    {
1055      if (segPhi)
1056      {
1057        // Use inner phi tolerant boundary -> if on tolerant
1058        // phi boundaries, phi intersect code handles leaving/entering checks
1059
1060        cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1061
1062        if (cosPsi>=cosHDPhiIT)
1063        { 
1064          // inside radii, delta r -ve, inside phi
1065
1066          if (segTheta)
1067          {
1068            if ( (pTheta >= tolSTheta + kAngTolerance)
1069              && (pTheta <= tolETheta - kAngTolerance) )
1070            {
1071              return snxt=0;
1072            }
1073          }
1074          else    // strictly inside Theta in both cases
1075          {
1076            return snxt=0;
1077          }
1078        }
1079      }
1080      else
1081      {
1082        if ( segTheta )
1083        {
1084          if ( (pTheta >= tolSTheta + kAngTolerance)
1085            && (pTheta <= tolETheta - kAngTolerance) )
1086          {
1087            return snxt=0;
1088          }
1089        }
1090        else   // strictly inside Theta in both cases
1091        {
1092          return snxt=0;
1093        }
1094      }
1095    }
1096  }
1097
1098  // Inner spherical shell intersection
1099  // - Always farthest root, because would have passed through outer
1100  //   surface first.
1101  // - Tolerant check for if travelling through solid
1102
1103  if (fRmin)
1104  {
1105    c  = rad2 - fRmin*fRmin ;
1106    d2 = pDotV3d*pDotV3d - c ;
1107
1108    // Within tolerance inner radius of inner G4Sphere
1109    // Check for immediate entry/already inside and travelling outwards
1110
1111    // if (c >- kRadTolerance*0.5 && pDotV3d >= 0 && rad2 < tolIRMin2 )
1112
1113    if ( c > -kRadTolerance*0.5 && rad2 < tolIRMin2 && 
1114         ( d2 < fRmin*kCarTolerance || pDotV3d >= 0 ) )
1115    {
1116      if (segPhi)
1117      {
1118        // Use inner phi tolerant boundary -> if on tolerant
1119        // phi boundaries, phi intersect code handles leaving/entering checks
1120
1121        cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
1122        if (cosPsi >= cosHDPhiIT)
1123        { 
1124          // inside radii, delta r -ve, inside phi
1125          //
1126          if (segTheta)
1127          {
1128            if ( (pTheta >= tolSTheta + kAngTolerance)
1129              && (pTheta <= tolETheta - kAngTolerance) )
1130            {
1131              return snxt=0;
1132            }
1133          }
1134          else
1135          {
1136            return snxt = 0 ;
1137          }
1138        }
1139      }
1140      else
1141      {
1142        if (segTheta)
1143        {
1144          if ( (pTheta >= tolSTheta + kAngTolerance)
1145            && (pTheta <= tolETheta - kAngTolerance) )
1146          {
1147            return snxt = 0 ;
1148          }
1149        }
1150        else
1151        {
1152          return snxt=0;
1153        }
1154      }
1155    }
1156    else   // Not special tolerant case
1157    {
1158      //  d2 = pDotV3d*pDotV3d - c ;
1159
1160      if (d2 >= 0)
1161      {
1162        s = -pDotV3d + std::sqrt(d2) ;
1163        if ( s >= kRadTolerance*0.5 )  // It was >= 0 ??
1164        {
1165          xi   = p.x() + s*v.x() ;
1166          yi   = p.y() + s*v.y() ;
1167          rhoi = std::sqrt(xi*xi+yi*yi) ;
1168
1169          if ( segPhi && rhoi )   // Check phi intersection
1170          {
1171            cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
1172
1173            if (cosPsi >= cosHDPhiOT)
1174            {
1175              if (segTheta)  // Check theta intersection
1176              {
1177                zi = p.z() + s*v.z() ;
1178
1179                // rhoi & zi can never both be 0
1180                // (=>intersect at origin =>fRmax=0)
1181                //
1182                iTheta = std::atan2(rhoi,zi) ;
1183                if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
1184                {
1185                  snxt = s ;
1186                }
1187              }
1188              else
1189              {
1190                snxt=s;
1191              }
1192            }
1193          }
1194          else
1195          {
1196            if (segTheta)   // Check theta intersection
1197            {
1198              zi = p.z() + s*v.z() ;
1199
1200              // rhoi & zi can never both be 0
1201              // (=>intersect at origin => fRmax=0 !)
1202              //
1203              iTheta = std::atan2(rhoi,zi) ;
1204              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1205              {
1206                snxt = s ;
1207              }
1208            }
1209            else
1210            {
1211              snxt=s;
1212            }
1213          }
1214        }
1215      }
1216    }
1217  }
1218
1219  // Phi segment intersection
1220  //
1221  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1222  //
1223  // o NOTE: Large duplication of code between sphi & ephi checks
1224  //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1225  //            intersection check <=0 -> >=0
1226  //         -> Should use some form of loop Construct
1227  //
1228  if ( segPhi )
1229  {
1230    // First phi surface (`S'tarting phi)
1231
1232    sinSPhi = std::sin(fSPhi) ;
1233    cosSPhi = std::cos(fSPhi) ;
1234
1235    // Comp = Component in outwards normal dirn
1236    //
1237    Comp    = v.x()*sinSPhi - v.y()*cosSPhi  ;
1238                   
1239    if ( Comp < 0 )
1240    {
1241      Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
1242
1243      if (Dist < kCarTolerance*0.5)
1244      {
1245        s = Dist/Comp ;
1246
1247        if (s < snxt)
1248        {
1249          if ( s > 0 )
1250          {
1251            xi    = p.x() + s*v.x() ;
1252            yi    = p.y() + s*v.y() ;
1253            zi    = p.z() + s*v.z() ;
1254            rhoi2 = xi*xi + yi*yi   ;
1255            radi2 = rhoi2 + zi*zi   ;
1256          }
1257          else
1258          {
1259            s     = 0     ;
1260            xi    = p.x() ;
1261            yi    = p.y() ;
1262            zi    = p.z() ;
1263            rhoi2 = rho2  ;
1264            radi2 = rad2  ;
1265          }
1266          if ( (radi2 <= tolORMax2)
1267            && (radi2 >= tolORMin2)
1268            && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
1269          {
1270            // Check theta intersection
1271            // rhoi & zi can never both be 0
1272            // (=>intersect at origin =>fRmax=0)
1273            //
1274            if ( segTheta )
1275            {
1276              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1277              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1278              {
1279                // r and theta intersections good
1280                // - check intersecting with correct half-plane
1281
1282                if ((yi*cosCPhi-xi*sinCPhi) <= 0)
1283                {
1284                  snxt = s ;
1285                }
1286              }
1287            }
1288            else
1289            {
1290              snxt = s ;
1291            }
1292          }
1293        }
1294      }
1295    }
1296
1297    // Second phi surface (`E'nding phi)
1298
1299    ePhi    = fSPhi + fDPhi ;
1300    sinEPhi = std::sin(ePhi)     ;
1301    cosEPhi = std::cos(ePhi)     ;
1302
1303    // Compnent in outwards normal dirn
1304
1305    Comp    = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
1306       
1307    if (Comp < 0)
1308    {
1309      Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
1310      if ( Dist < kCarTolerance*0.5 )
1311      {
1312        s = Dist/Comp ;
1313
1314        if ( s < snxt )
1315        {
1316          if (s > 0)
1317          {
1318            xi    = p.x() + s*v.x() ;
1319            yi    = p.y() + s*v.y() ;
1320            zi    = p.z() + s*v.z() ;
1321            rhoi2 = xi*xi + yi*yi   ;
1322            radi2 = rhoi2 + zi*zi   ;
1323          }
1324          else
1325          {
1326            s     = 0     ;
1327            xi    = p.x() ;
1328            yi    = p.y() ;
1329            zi    = p.z() ;
1330            rhoi2 = rho2  ;
1331            radi2 = rad2  ;
1332          } if ( (radi2 <= tolORMax2)
1333            && (radi2 >= tolORMin2)
1334            && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
1335          {
1336            // Check theta intersection
1337            // rhoi & zi can never both be 0
1338            // (=>intersect at origin =>fRmax=0)
1339            //
1340            if ( segTheta )
1341            {
1342              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1343              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1344              {
1345                // r and theta intersections good
1346                // - check intersecting with correct half-plane
1347
1348                if ((yi*cosCPhi-xi*sinCPhi) >= 0)
1349                {
1350                  snxt = s ;
1351                }
1352              }
1353            }
1354            else
1355            {
1356              snxt = s ;
1357            }
1358          }
1359        }
1360      }
1361    }
1362  }
1363
1364  // Theta segment intersection
1365
1366  if ( segTheta )
1367  {
1368
1369    // Intersection with theta surfaces
1370    // Known failure cases:
1371    // o  Inside tolerance of stheta surface, skim
1372    //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
1373    //
1374    //    To solve: Check 2nd root of etheta surface in addition to stheta
1375    //
1376    // o  start/end theta is exactly pi/2
1377    // Intersections with cones
1378    //
1379    // Cone equation: x^2+y^2=z^2tan^2(t)
1380    //
1381    // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1382    //
1383    // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t))
1384    //       + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1385    //
1386    // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1387
1388    tanSTheta  = std::tan(fSTheta)         ;
1389    tanSTheta2 = tanSTheta*tanSTheta  ;
1390    tanETheta  = std::tan(fSTheta+fDTheta) ;
1391    tanETheta2 = tanETheta*tanETheta  ;
1392     
1393    if (fSTheta)
1394    {
1395      dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
1396    }
1397    else
1398    {
1399      dist2STheta = kInfinity ;
1400    }
1401    if ( fSTheta + fDTheta < pi )
1402    {
1403      dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
1404    }
1405      else
1406    {
1407      dist2ETheta=kInfinity;
1408    }     
1409    if ( pTheta < tolSTheta) // dist2STheta<-kRadTolerance*0.5 && dist2ETheta>0)
1410    {
1411      // Inside (theta<stheta-tol) s theta cone
1412      // First root of stheta cone, second if first root -ve
1413
1414      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1415      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1416       
1417      b  = t2/t1 ;
1418      c  = dist2STheta/t1 ;
1419      d2 = b*b - c ;
1420
1421      if ( d2 >= 0 )
1422      {
1423        d = std::sqrt(d2) ;
1424        s = -b - d ;    // First root
1425
1426        if ( s < 0 )
1427        {
1428          s=-b+d;    // Second root
1429        }
1430        if (s >= 0 && s < snxt)
1431        {
1432          xi    = p.x() + s*v.x() ;
1433          yi    = p.y() + s*v.y() ;
1434          zi    = p.z() + s*v.z() ;
1435          rhoi2 = xi*xi + yi*yi   ;
1436          radi2 = rhoi2 + zi*zi   ;
1437          if ( (radi2 <= tolORMax2)
1438            && (radi2 >= tolORMin2)
1439            && (zi*(fSTheta - halfpi) <= 0) )
1440          {
1441            if ( segPhi && rhoi2 )  // Check phi intersection
1442            {
1443              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1444              if (cosPsi >= cosHDPhiOT)
1445              {
1446                snxt = s ;
1447              }
1448            }
1449            else
1450            {
1451              snxt = s ;
1452            }
1453          }
1454        }
1455      }
1456
1457      // Possible intersection with ETheta cone.
1458      // Second >= 0 root should be considered
1459       
1460      if ( fSTheta + fDTheta < pi )
1461      {
1462        t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1463        t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1464       
1465        b  = t2/t1 ;
1466        c  = dist2ETheta/t1 ;
1467        d2 = b*b - c ;
1468
1469        if (d2 >= 0)
1470        {
1471          d = std::sqrt(d2) ;
1472          s = -b + d ;    // Second root
1473
1474          if (s >= 0 && s < snxt)
1475          {
1476            xi    = p.x() + s*v.x() ;
1477            yi    = p.y() + s*v.y() ;
1478            zi    = p.z() + s*v.z() ;
1479            rhoi2 = xi*xi + yi*yi   ;
1480            radi2 = rhoi2 + zi*zi   ;
1481
1482            if ( (radi2 <= tolORMax2)
1483              && (radi2 >= tolORMin2)
1484              && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
1485            {
1486              if (segPhi && rhoi2)   // Check phi intersection
1487              {
1488                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1489                if (cosPsi >= cosHDPhiOT)
1490                {
1491                  snxt = s ;
1492                }
1493              }
1494              else
1495              {
1496                snxt = s ;
1497              }
1498            }
1499          }
1500        }
1501      }
1502    } 
1503    else if (pTheta > tolETheta) 
1504    { // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
1505      // Inside (theta>etheta+tol) e theta cone
1506      // First root of etheta cone, second if first root `imaginary'
1507
1508      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1509      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1510       
1511      b  = t2/t1 ;
1512      c  = dist2ETheta/t1 ;
1513      d2 = b*b - c ;
1514
1515      if (d2 >= 0)
1516      {
1517        d = std::sqrt(d2) ;
1518        s = -b - d ;    // First root
1519        if (s < 0)
1520        {
1521          s = -b + d ;           // second root
1522        }
1523        if (s >= 0 && s < snxt)
1524        {
1525          xi    = p.x() + s*v.x() ;
1526          yi    = p.y() + s*v.y() ;
1527          zi    = p.z() + s*v.z() ;
1528          rhoi2 = xi*xi + yi*yi   ;
1529          radi2 = rhoi2 + zi*zi   ;
1530
1531          if ( (radi2 <= tolORMax2)
1532            && (radi2 >= tolORMin2) 
1533            && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
1534          {
1535            if (segPhi && rhoi2)  // Check phi intersection
1536            {
1537              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1538              if (cosPsi >= cosHDPhiOT)
1539              {
1540                snxt = s ;
1541              }
1542            }
1543            else
1544            {
1545              snxt = s ;
1546            }
1547          }
1548        }
1549      }
1550
1551      // Possible intersection with STheta cone.
1552      // Second >= 0 root should be considered
1553       
1554      if ( fSTheta )
1555      {
1556        t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1557        t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1558
1559        b  = t2/t1 ;
1560        c  = dist2STheta/t1 ;
1561        d2 = b*b - c ;
1562
1563        if (d2 >= 0)
1564        {
1565          d = std::sqrt(d2) ;
1566          s = -b + d ;    // Second root
1567
1568          if ( (s >= 0) && (s < snxt) )
1569          {
1570            xi    = p.x() + s*v.x() ;
1571            yi    = p.y() + s*v.y() ;
1572            zi    = p.z() + s*v.z() ;
1573            rhoi2 = xi*xi + yi*yi   ;
1574            radi2 = rhoi2 + zi*zi   ;
1575
1576            if ( (radi2 <= tolORMax2)
1577              && (radi2 >= tolORMin2)
1578              && (zi*(fSTheta - halfpi) <= 0) )
1579            {
1580              if (segPhi && rhoi2)   // Check phi intersection
1581              {
1582                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1583                if (cosPsi >= cosHDPhiOT)
1584                {
1585                  snxt = s ;
1586                }
1587              }
1588              else
1589              {
1590                snxt = s ;
1591              }
1592            }
1593          }
1594        }
1595      } 
1596    }     
1597    else if ( (pTheta <tolSTheta + kAngTolerance)
1598           && (fSTheta > kAngTolerance) )
1599    {
1600      // In tolerance of stheta
1601      // If entering through solid [r,phi] => 0 to in
1602      // else try 2nd root
1603
1604      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1605      if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<pi*.5)
1606        || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>pi*.5)
1607        || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==pi*.5) )
1608      {
1609        if (segPhi && rho2)  // Check phi intersection
1610        {
1611          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1612          if (cosPsi >= cosHDPhiIT)
1613          {
1614            return 0 ;
1615          }
1616        }
1617        else
1618        {
1619          return 0 ;
1620        }
1621      }
1622
1623      // Not entering immediately/travelling through
1624
1625      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1626      b  = t2/t1 ;
1627      c  = dist2STheta/t1 ;
1628      d2 = b*b - c ;
1629
1630      if (d2 >= 0)
1631      {
1632        d = std::sqrt(d2) ;
1633        s = -b + d ;
1634        if ( (s >= kCarTolerance*0.5) && (s < snxt) && (fSTheta < pi*0.5) )
1635        {
1636          xi    = p.x() + s*v.x() ;
1637          yi    = p.y() + s*v.y() ;
1638          zi    = p.z() + s*v.z() ;
1639          rhoi2 = xi*xi + yi*yi   ;
1640          radi2 = rhoi2 + zi*zi   ;
1641
1642          if ( (radi2 <= tolORMax2)
1643            && (radi2 >= tolORMin2)
1644            && (zi*(fSTheta - halfpi) <= 0) )
1645          {
1646            if ( segPhi && rhoi2 )    // Check phi intersection
1647            {
1648              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1649              if ( cosPsi >= cosHDPhiOT )
1650              {
1651                snxt = s ;
1652              }
1653            }
1654            else
1655            {
1656              snxt = s ;
1657            }
1658          }
1659        }
1660      }
1661    }   
1662    else if ( (pTheta > tolETheta - kAngTolerance)
1663           && ((fSTheta + fDTheta) < pi-kAngTolerance) )   
1664    {
1665
1666      // In tolerance of etheta
1667      // If entering through solid [r,phi] => 0 to in
1668      // else try 2nd root
1669
1670      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1671
1672      if (
1673    (t2<0    && (fSTheta+fDTheta) <pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
1674 || (t2>=0   && (fSTheta+fDTheta) >pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
1675 || (v.z()>0 && (fSTheta+fDTheta)==pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
1676         )
1677      {
1678        if (segPhi && rho2)   // Check phi intersection
1679        {
1680          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1681          if (cosPsi >= cosHDPhiIT)
1682          {
1683            return 0 ;
1684          }
1685        }
1686        else
1687        {
1688          return 0 ;
1689        }
1690      }
1691
1692      // Not entering immediately/travelling through
1693
1694      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1695      b  = t2/t1 ;
1696      c  = dist2ETheta/t1 ;
1697      d2 = b*b - c ;
1698
1699      if (d2 >= 0)
1700      {
1701        d = std::sqrt(d2) ;
1702        s = -b + d ;
1703       
1704        if ( (s >= kCarTolerance*0.5)
1705          && (s < snxt) && ((fSTheta + fDTheta) > pi*0.5) )
1706        {
1707          xi    = p.x() + s*v.x() ;
1708          yi    = p.y() + s*v.y() ;
1709          zi    = p.z() + s*v.z() ;
1710          rhoi2 = xi*xi + yi*yi   ;
1711          radi2 = rhoi2 + zi*zi   ;
1712
1713          if ( (radi2 <= tolORMax2)
1714            && (radi2 >= tolORMin2)
1715            && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
1716          {
1717            if (segPhi && rhoi2)   // Check phi intersection
1718            {
1719              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1720              if (cosPsi>=cosHDPhiOT)
1721              {
1722                snxt = s ;
1723              }
1724            }
1725            else
1726            {
1727              snxt = s ;
1728            }
1729          }
1730        }
1731      }   
1732    } 
1733    else
1734    {
1735      // stheta+tol<theta<etheta-tol
1736      // For BOTH stheta & etheta check 2nd root for validity [r,phi]
1737
1738      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1739      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1740
1741      b  = t2/t1;
1742      c  = dist2STheta/t1 ;
1743      d2 = b*b - c ;
1744
1745      if (d2 >= 0)
1746      {
1747        d = std::sqrt(d2) ;
1748        s = -b + d ;    // second root
1749
1750        if (s >= 0 && s < snxt)
1751        {
1752          xi    = p.x() + s*v.x() ;
1753          yi    = p.y() + s*v.y() ;
1754          zi    = p.z() + s*v.z() ;
1755          rhoi2 = xi*xi + yi*yi   ;
1756          radi2 = rhoi2 + zi*zi   ;
1757
1758          if ( (radi2 <= tolORMax2)
1759            && (radi2 >= tolORMin2)
1760            && (zi*(fSTheta - halfpi) <= 0) )
1761          {
1762            if (segPhi && rhoi2)   // Check phi intersection
1763            {
1764              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1765              if (cosPsi >= cosHDPhiOT)
1766              {
1767                snxt = s ;
1768              }
1769            }
1770            else
1771            {
1772              snxt = s ;
1773            }
1774          }
1775        }
1776      }       
1777      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1778      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1779       
1780      b  = t2/t1 ;
1781      c  = dist2ETheta/t1 ;
1782      d2 = b*b - c ;
1783
1784      if (d2 >= 0)
1785      {
1786        d = std::sqrt(d2) ;
1787        s = -b + d;    // second root
1788
1789        if (s >= 0 && s < snxt)
1790        {
1791          xi    = p.x() + s*v.x() ;
1792          yi    = p.y() + s*v.y() ;
1793          zi    = p.z() + s*v.z() ;
1794          rhoi2 = xi*xi + yi*yi   ;
1795          radi2 = rhoi2 + zi*zi   ;
1796
1797          if ( (radi2 <= tolORMax2)
1798            && (radi2 >= tolORMin2)
1799            && (zi*(fSTheta + fDTheta - halfpi) <= 0) )
1800          {
1801            if (segPhi && rhoi2)   // Check phi intersection
1802            {
1803              cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1804              if ( cosPsi >= cosHDPhiOT )
1805              {
1806                snxt=s;
1807              }
1808            }
1809            else
1810            {
1811              snxt = s ;
1812            }
1813          }
1814        }
1815      }
1816    } 
1817  }
1818  return snxt;
1819}
1820
1821//////////////////////////////////////////////////////////////////////
1822//
1823// Calculate distance (<= actual) to closest surface of shape from outside
1824// - Calculate distance to radial planes
1825// - Only to phi planes if outside phi extent
1826// - Only to theta planes if outside theta extent
1827// - Return 0 if point inside
1828
1829G4double G4Sphere::DistanceToIn( const G4ThreeVector& p ) const
1830{
1831  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
1832  G4double rho2,rad,rho;
1833  G4double phiC,cosPhiC,sinPhiC,cosPsi,ePhi;
1834  G4double pTheta,dTheta1,dTheta2;
1835  rho2=p.x()*p.x()+p.y()*p.y();
1836  rad=std::sqrt(rho2+p.z()*p.z());
1837  rho=std::sqrt(rho2);
1838
1839  //
1840  // Distance to r shells
1841  //   
1842  if (fRmin)
1843  {
1844    safeRMin=fRmin-rad;
1845    safeRMax=rad-fRmax;
1846    if (safeRMin>safeRMax)
1847    {
1848      safe=safeRMin;
1849    }
1850    else
1851    {
1852      safe=safeRMax;
1853    }
1854  }
1855  else
1856  {
1857    safe=rad-fRmax;
1858  }
1859
1860  //
1861  // Distance to phi extent
1862  //
1863  if (fDPhi<twopi&&rho)
1864  {
1865    phiC=fSPhi+fDPhi*0.5;
1866    cosPhiC=std::cos(phiC);
1867    sinPhiC=std::sin(phiC);
1868
1869    // Psi=angle from central phi to point
1870    //
1871    cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho;
1872    if (cosPsi<std::cos(fDPhi*0.5))
1873    {
1874      // Point lies outside phi range
1875      //
1876      if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1877      {
1878        safePhi=std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1879      }
1880      else
1881      {
1882        ePhi=fSPhi+fDPhi;
1883        safePhi=std::fabs(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
1884      }
1885      if (safePhi>safe) safe=safePhi;
1886    }
1887  }
1888  //
1889  // Distance to Theta extent
1890  //   
1891  if ((rad!=0.0) && (fDTheta<pi))
1892  {
1893    pTheta=std::acos(p.z()/rad);
1894    if (pTheta<0) pTheta+=pi;
1895    dTheta1=fSTheta-pTheta;
1896    dTheta2=pTheta-(fSTheta+fDTheta);
1897    if (dTheta1>dTheta2)
1898    {
1899      if (dTheta1>=0)             // WHY ???????????
1900      {
1901        safeTheta=rad*std::sin(dTheta1);
1902        if (safe<=safeTheta)
1903        {
1904          safe=safeTheta;
1905        }
1906      }
1907    }
1908    else
1909    {
1910      if (dTheta2>=0)
1911      {
1912        safeTheta=rad*std::sin(dTheta2);
1913        if (safe<=safeTheta)
1914        {
1915          safe=safeTheta;
1916        }
1917      }
1918    }
1919  }
1920
1921  if (safe<0) safe=0;
1922  return safe;
1923}
1924
1925/////////////////////////////////////////////////////////////////////
1926//
1927// Calculate distance to surface of shape from `inside', allowing for tolerance
1928// - Only Calc rmax intersection if no valid rmin intersection
1929
1930G4double G4Sphere::DistanceToOut( const G4ThreeVector& p,
1931                                  const G4ThreeVector& v,
1932                                  const G4bool calcNorm,
1933                                        G4bool *validNorm,
1934                                        G4ThreeVector *n   ) const
1935{
1936  G4double snxt = kInfinity;     // snxt is default return value
1937  G4double sphi= kInfinity,stheta= kInfinity;
1938  ESide side=kNull,sidephi=kNull,sidetheta=kNull; 
1939
1940  G4double t1,t2;
1941  G4double b,c,d;
1942
1943  // Variables for phi intersection:
1944
1945  G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi;
1946  G4double cPhi,sinCPhi,cosCPhi;
1947  G4double pDistS,compS,pDistE,compE,sphi2,vphi;
1948   
1949  G4double rho2,rad2,pDotV2d,pDotV3d,pTheta;
1950
1951  G4double tolSTheta=0.,tolETheta=0.;
1952  G4double xi,yi,zi;      // Intersection point
1953
1954  // G4double Comp; // Phi intersection
1955
1956  G4bool segPhi;        // Phi flag and precalcs
1957  G4double hDPhi,hDPhiOT,hDPhiIT; 
1958  G4double cosHDPhiOT,cosHDPhiIT;
1959
1960  G4bool segTheta;                             // Theta flag and precals
1961  G4double tanSTheta=0.,tanETheta, rhoSecTheta;
1962  G4double tanSTheta2=0.,tanETheta2=0.;
1963  G4double dist2STheta,dist2ETheta;
1964  G4double d2,s;
1965
1966  // General Precalcs
1967
1968  rho2=p.x()*p.x()+p.y()*p.y();
1969  rad2=rho2+p.z()*p.z();
1970  //  G4double rad=std::sqrt(rad2);
1971
1972  pTheta=std::atan2(std::sqrt(rho2),p.z());
1973
1974  pDotV2d=p.x()*v.x()+p.y()*v.y();
1975  pDotV3d=pDotV2d+p.z()*v.z();
1976
1977  // Set phi divided flag and precalcs
1978
1979  if(fDPhi<twopi)
1980  {
1981    segPhi=true;
1982    hDPhi=0.5*fDPhi;    // half delta phi
1983    cPhi=fSPhi+hDPhi;;
1984    hDPhiOT=hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi
1985    hDPhiIT=hDPhi-0.5*kAngTolerance;
1986    sinCPhi=std::sin(cPhi);
1987    cosCPhi=std::cos(cPhi);
1988    cosHDPhiOT=std::cos(hDPhiOT);
1989    cosHDPhiIT=std::cos(hDPhiIT);
1990  }
1991  else
1992  {
1993    segPhi=false;
1994  }
1995
1996  // Theta precalcs
1997   
1998  if (fDTheta < pi)
1999  {
2000    segTheta=true;
2001    tolSTheta=fSTheta-kAngTolerance*0.5;
2002    tolETheta=fSTheta+fDTheta+kAngTolerance*0.5;
2003  }
2004  else
2005  {
2006    segTheta=false;
2007  }
2008   
2009  // Radial Intersections from G4Sphere::DistanceToIn
2010  //
2011  // Outer spherical shell intersection
2012  // - Only if outside tolerant fRmax
2013  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
2014  // - No intersect -> no intersection with G4Sphere
2015  //
2016  // Shell eqn: x^2+y^2+z^2=RSPH^2
2017  //
2018  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
2019  //
2020  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
2021  // =>      rad2        +2s(pDotV3d)       +s^2                =R^2
2022  //
2023  // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
2024  //
2025  // const G4double  fractionTolerance = 1.0e-12;
2026  const G4double  flexRadMaxTolerance = // kRadTolerance;
2027    std::max(kRadTolerance, fEpsilon * fRmax);
2028
2029  const G4double  Rmax_plus = fRmax + flexRadMaxTolerance*0.5;
2030  const G4double  flexRadMinTolerance = std::max(kRadTolerance, 
2031                     fEpsilon * fRmin);
2032  const G4double  Rmin_minus= (fRmin > 0) ? fRmin-flexRadMinTolerance*0.5 : 0 ;
2033
2034  if(rad2 <= Rmax_plus*Rmax_plus && rad2 >= Rmin_minus*Rmin_minus)
2035    //  if(rad <= Rmax_plus && rad >= Rmin_minus)
2036  {
2037    c = rad2 - fRmax*fRmax;
2038
2039    if (c < flexRadMaxTolerance*fRmax) 
2040    {
2041      // Within tolerant Outer radius
2042      //
2043      // The test is
2044      //     rad  - fRmax < 0.5*kRadTolerance
2045      // =>  rad  < fRmax + 0.5*kRadTol
2046      // =>  rad2 < (fRmax + 0.5*kRadTol)^2
2047      // =>  rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
2048      // =>  rad2 - fRmax^2    <~    fRmax*kRadTol
2049
2050      d2 = pDotV3d*pDotV3d - c;
2051
2052      if( (c >- flexRadMaxTolerance*fRmax)       // on tolerant surface
2053       && ((pDotV3d >=0) || (d2 < 0)) )          // leaving outside from Rmax
2054                                                 // not re-entering
2055      {
2056        if(calcNorm)
2057        {
2058          *validNorm = true ;
2059          *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
2060        }
2061        return snxt = 0;
2062      }
2063      else 
2064      {
2065        snxt=-pDotV3d+std::sqrt(d2);    // second root since inside Rmax
2066        side = kRMax ; 
2067      }
2068    }
2069
2070    // Inner spherical shell intersection:
2071    // Always first >=0 root, because would have passed
2072    // from outside of Rmin surface .
2073
2074    if (fRmin)
2075    {
2076      c  = rad2 - fRmin*fRmin;
2077      d2 = pDotV3d*pDotV3d - c;
2078
2079      if (c >- flexRadMinTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
2080      {
2081        if( c < flexRadMinTolerance*fRmin && 
2082            d2 >= flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin
2083        {
2084          if(calcNorm)
2085          {
2086            *validNorm = false ;   // Rmin surface is concave
2087          }
2088          return snxt = 0 ;
2089        }
2090        else
2091        { 
2092          if (d2 >= 0)
2093          {
2094            s = -pDotV3d-std::sqrt(d2) ;
2095            if (s>=0)     // Always intersect Rmin first
2096            {
2097              snxt = s ;
2098              side = kRMin ;
2099            }
2100          }
2101        }
2102      }
2103    }
2104  }
2105
2106  // Theta segment intersection
2107
2108  if (segTheta)
2109  {
2110    // Intersection with theta surfaces
2111    //
2112    // Known failure cases:
2113    // o  Inside tolerance of stheta surface, skim
2114    //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
2115    //
2116    //    To solve: Check 2nd root of etheta surface in addition to stheta
2117    //
2118    // o  start/end theta is exactly pi/2
2119    //
2120    // Intersections with cones
2121    //
2122    // Cone equation: x^2+y^2=z^2tan^2(t)
2123    //
2124    // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
2125    //
2126    // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t))
2127    //       + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0
2128    //
2129    // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
2130    //
2131    tanSTheta=std::tan(fSTheta);
2132    tanSTheta2=tanSTheta*tanSTheta;
2133    tanETheta=std::tan(fSTheta+fDTheta);
2134    tanETheta2=tanETheta*tanETheta;
2135     
2136    if (fSTheta)
2137    {
2138      dist2STheta=rho2-p.z()*p.z()*tanSTheta2;
2139    }
2140    else
2141    {
2142      dist2STheta = kInfinity;
2143    }
2144    if (fSTheta + fDTheta < pi)
2145    {
2146      dist2ETheta = rho2-p.z()*p.z()*tanETheta2;
2147    }
2148    else
2149    {
2150      dist2ETheta = kInfinity ;
2151    }
2152    if (pTheta > tolSTheta && pTheta < tolETheta)   // Inside theta 
2153    {
2154      // In tolerance of STheta and possible leaving out to small thetas N-
2155
2156      if(pTheta < tolSTheta + kAngTolerance  && fSTheta > kAngTolerance) 
2157      {
2158        t2=pDotV2d-p.z()*v.z()*tanSTheta2 ; // =(VdotN+)*rhoSecSTheta
2159
2160        if( fSTheta < pi*0.5 && t2 < 0)
2161        {
2162          if(calcNorm) *validNorm = false ;
2163          return snxt = 0 ;
2164        }
2165        else if(fSTheta > pi*0.5 && t2 >= 0)
2166        {
2167          if(calcNorm)
2168          {
2169            rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2)) ;
2170            *validNorm = true ;
2171            *n = G4ThreeVector(-p.x()/rhoSecTheta,   // N-
2172                               -p.y()/rhoSecTheta,
2173                               tanSTheta/std::sqrt(1+tanSTheta2) ) ;
2174          }
2175          return snxt = 0 ;
2176        }
2177        else if( fSTheta == pi*0.5 && v.z() > 0)
2178        {
2179          if(calcNorm)
2180          {
2181            *validNorm = true ;
2182            *n = G4ThreeVector(0,0,1) ;
2183          }
2184          return snxt = 0 ;
2185        }
2186      }
2187
2188      // In tolerance of ETheta and possible leaving out to larger thetas N+
2189
2190      if ( (pTheta  > tolETheta - kAngTolerance)
2191        && (( fSTheta + fDTheta) < pi - kAngTolerance) ) 
2192      {
2193        t2=pDotV2d-p.z()*v.z()*tanETheta2 ;
2194        if((fSTheta+fDTheta)>pi*0.5 && t2<0)
2195        {
2196          if(calcNorm) *validNorm = false ;
2197          return snxt = 0 ;
2198        }
2199        else if( (fSTheta+fDTheta) < pi*0.5 && t2 >= 0 )
2200        {
2201          if(calcNorm)
2202          {
2203            rhoSecTheta = std::sqrt(rho2*(1+tanETheta2)) ;
2204            *validNorm = true ;
2205            *n = G4ThreeVector( p.x()/rhoSecTheta,  // N+
2206                                p.y()/rhoSecTheta,
2207                                -tanETheta/std::sqrt(1+tanETheta2)  ) ; 
2208          }
2209          return snxt = 0 ;
2210        }
2211        else if( ( fSTheta+fDTheta) == pi*0.5 && v.z() < 0 )
2212        {
2213          if(calcNorm)
2214          {
2215            *validNorm = true ;
2216            *n = G4ThreeVector(0,0,-1) ;
2217          }
2218          return snxt = 0 ;
2219        }
2220      }
2221      if( fSTheta > 0 )
2222      {       
2223        // First root of fSTheta cone, second if first root -ve
2224
2225        t1 = 1-v.z()*v.z()*(1+tanSTheta2);
2226        t2 = pDotV2d-p.z()*v.z()*tanSTheta2;
2227       
2228        b  = t2/t1;
2229        c  = dist2STheta/t1;
2230        d2 = b*b - c ;
2231
2232        if ( d2 >= 0 )
2233        {
2234          d = std::sqrt(d2) ;
2235          s = -b - d ;    // First root
2236
2237          if ( s < 0 )
2238          {
2239            s = -b + d ;    // Second root
2240          }
2241          if (s > flexRadMaxTolerance*0.5 )   // && s<sr)
2242          {
2243            // check against double cone solution
2244            zi=p.z()+s*v.z();
2245            if (fSTheta<pi*0.5 && zi<0)
2246            {
2247              s = kInfinity ;  // wrong cone
2248            }
2249            if (fSTheta>pi*0.5 && zi>0)
2250            {
2251              s = kInfinity ;  // wrong cone
2252            }
2253            stheta = s ;
2254            sidetheta = kSTheta ;
2255          }
2256        }
2257      }
2258
2259      // Possible intersection with ETheta cone 
2260     
2261      if (fSTheta + fDTheta < pi)
2262      {
2263        t1 = 1-v.z()*v.z()*(1+tanETheta2);
2264        t2 = pDotV2d-p.z()*v.z()*tanETheta2;       
2265        b  = t2/t1;
2266        c  = dist2ETheta/t1;
2267        d2 = b*b-c ;
2268
2269        if ( d2 >= 0 )
2270        {
2271          d = std::sqrt(d2);
2272          s = -b - d ;          // First root
2273
2274          if ( s < 0 )
2275          {
2276            s=-b+d;    // Second root
2277          }
2278          if (s > flexRadMaxTolerance*0.5 && s < stheta )
2279          {
2280            // check against double cone solution
2281            zi=p.z()+s*v.z();
2282            if (fSTheta+fDTheta<pi*0.5 && zi<0)
2283            {
2284              s = kInfinity ;  // wrong cone
2285            }
2286            if (fSTheta+fDTheta>pi*0.5 && zi>0)
2287            {
2288              s = kInfinity ;  // wrong cone
2289            }
2290            if (s < stheta)
2291            {
2292              stheta = s ;
2293              sidetheta = kETheta ;
2294            }
2295          }
2296        }
2297      }
2298    } 
2299  }
2300
2301  // Phi Intersection
2302   
2303  if ( fDPhi < twopi)
2304  {
2305    sinSPhi=std::sin(fSPhi);
2306    cosSPhi=std::cos(fSPhi);
2307    ePhi=fSPhi+fDPhi;
2308    sinEPhi=std::sin(ePhi);
2309    cosEPhi=std::cos(ePhi);
2310    cPhi=fSPhi+fDPhi*0.5;
2311    sinCPhi=std::sin(cPhi);
2312    cosCPhi=std::cos(cPhi);
2313
2314    if ( p.x()||p.y() ) // Check if on z axis (rho not needed later)
2315    {
2316      // pDist -ve when inside
2317
2318      pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
2319      pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
2320
2321      // Comp -ve when in direction of outwards normal
2322
2323      compS   = -sinSPhi*v.x()+cosSPhi*v.y() ;
2324      compE   =  sinEPhi*v.x()-cosEPhi*v.y() ;
2325      sidephi = kNull ;
2326
2327      if ( pDistS <= 0 && pDistE <= 0 )
2328      {
2329        // Inside both phi *full* planes
2330
2331        if ( compS < 0 )
2332        {
2333          sphi = pDistS/compS ;
2334          xi   = p.x()+sphi*v.x() ;
2335          yi   = p.y()+sphi*v.y() ;
2336
2337          // Check intersecting with correct half-plane
2338          // (if not -> no intersect)
2339
2340          if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2341          {
2342            sphi=kInfinity;
2343          }
2344          else
2345          {
2346            sidephi = kSPhi ;
2347            if ( pDistS > -0.5*kCarTolerance) sphi =0 ; // Leave by sphi
2348          }
2349        }
2350        else sphi = kInfinity ;
2351
2352        if ( compE < 0 )
2353        {
2354          sphi2=pDistE/compE ;
2355          if (sphi2 < sphi) // Only check further if < starting phi intersection
2356          {
2357            xi = p.x()+sphi2*v.x() ;
2358            yi = p.y()+sphi2*v.y() ;
2359
2360            // Check intersecting with correct half-plane
2361 
2362            if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
2363            {
2364              sidephi = kEPhi ;
2365              if ( pDistE <= -0.5*kCarTolerance )
2366              {
2367                sphi=sphi2;
2368              }
2369              else 
2370              {
2371                sphi = 0 ;
2372              }
2373            }
2374          }
2375        }       
2376      }
2377      else if ( pDistS >= 0 && pDistE >= 0 ) // Outside both *full* phi planes
2378      {
2379        if ( pDistS <= pDistE )
2380        {
2381          sidephi = kSPhi ;
2382        }
2383        else
2384        {
2385          sidephi = kEPhi ;
2386        }
2387        if ( fDPhi > pi )
2388        {
2389          if ( compS < 0 && compE < 0 ) sphi = 0 ;
2390          else                          sphi = kInfinity ;
2391        }
2392        else
2393        {
2394          // if towards both >=0 then once inside (after error)
2395          // will remain inside
2396
2397          if ( compS >= 0 && compE >= 0 )
2398          {
2399            sphi=kInfinity;
2400          }
2401          else
2402          {
2403            sphi=0;
2404          }
2405        }   
2406      }
2407      else if ( pDistS > 0 && pDistE < 0 )
2408      {
2409        // Outside full starting plane, inside full ending plane
2410
2411        if ( fDPhi > pi )
2412        {
2413          if ( compE < 0 )
2414          {
2415            sphi = pDistE/compE ;
2416            xi   = p.x() + sphi*v.x() ;
2417            yi   = p.y() + sphi*v.y() ;
2418
2419            // Check intersection in correct half-plane
2420            // (if not -> not leaving phi extent)
2421            //
2422            if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
2423            {
2424              sphi = kInfinity ;
2425            }
2426            else // Leaving via Ending phi
2427            {
2428              sidephi = kEPhi ;
2429              if ( pDistE > -0.5*kCarTolerance ) sphi = 0. ;
2430            }
2431          }
2432          else
2433          {
2434            sphi = kInfinity ;
2435          }
2436        }
2437        else
2438        {
2439          if ( compS >= 0 )
2440          {
2441            if ( compE < 0 )
2442            {           
2443              sphi = pDistE/compE ;
2444              xi   = p.x() + sphi*v.x() ;
2445              yi   = p.y() + sphi*v.y() ;
2446
2447              // Check intersection in correct half-plane
2448              // (if not -> remain in extent)
2449              //
2450              if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
2451              {
2452                sphi=kInfinity;
2453              }
2454              else // otherwise leaving via Ending phi
2455              {
2456                sidephi = kEPhi ;
2457              }
2458            }
2459            else sphi=kInfinity;
2460          }
2461          else // leaving immediately by starting phi
2462          {
2463            sidephi = kSPhi ;
2464            sphi    = 0 ;
2465          }
2466        }
2467      }
2468      else
2469      {
2470        // Must be pDistS < 0 && pDistE > 0
2471        // Inside full starting plane, outside full ending plane
2472
2473        if ( fDPhi > pi )
2474        {
2475          if ( compS < 0 )
2476          {
2477            sphi=pDistS/compS;
2478            xi=p.x()+sphi*v.x();
2479            yi=p.y()+sphi*v.y();
2480
2481            // Check intersection in correct half-plane
2482            // (if not -> not leaving phi extent)
2483            //
2484            if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2485            {
2486              sphi = kInfinity ;
2487            }
2488           else  // Leaving via Starting phi
2489           {
2490              sidephi = kSPhi ;
2491             if ( pDistS > -0.5*kCarTolerance ) sphi = 0 ;
2492            }
2493          }
2494          else
2495          {
2496            sphi = kInfinity ;
2497          }
2498        }
2499        else
2500        {
2501          if ( compE >= 0 )
2502          {
2503            if ( compS < 0 )
2504            {
2505              sphi = pDistS/compS ;
2506              xi   = p.x()+sphi*v.x() ;
2507              yi   = p.y()+sphi*v.y() ;
2508
2509              // Check intersection in correct half-plane
2510              // (if not -> remain in extent)
2511              //
2512              if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2513              {
2514                sphi = kInfinity ;
2515              }
2516              else // otherwise leaving via Starting phi
2517              {
2518                sidephi = kSPhi ;
2519              }
2520            }
2521            else
2522            {
2523              sphi = kInfinity ;
2524            }
2525          }
2526          else // leaving immediately by ending
2527          {
2528            sidephi = kEPhi ;
2529            sphi    = 0     ;
2530          }
2531        }
2532      }     
2533    }
2534    else
2535    {
2536      // On z axis + travel not || to z axis -> if phi of vector direction
2537      // within phi of shape, Step limited by rmax, else Step =0
2538
2539      if ( v.x() || v.y() )
2540      {
2541        vphi = std::atan2(v.y(),v.x()) ;
2542        if ( fSPhi < vphi && vphi < fSPhi + fDPhi )
2543        {
2544          sphi=kInfinity;
2545        }
2546        else
2547        {
2548          sidephi = kSPhi ; // arbitrary
2549          sphi    = 0     ;
2550        }
2551      }
2552      else  // travel along z - no phi intersaction
2553      {
2554        sphi = kInfinity ;
2555      }
2556    }
2557    if ( sphi < snxt )  // Order intersecttions
2558    {
2559      snxt = sphi ;
2560      side = sidephi ;
2561    }
2562  }
2563  if (stheta < snxt ) // Order intersections
2564  {
2565    snxt = stheta ;
2566    side = sidetheta ;
2567  }
2568
2569  if (calcNorm)    // Output switch operator
2570  {
2571    switch( side )
2572    {
2573      case kRMax:
2574        xi=p.x()+snxt*v.x();
2575        yi=p.y()+snxt*v.y();
2576        zi=p.z()+snxt*v.z();
2577        *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
2578        *validNorm=true;
2579        break;
2580      case kRMin:
2581        *validNorm=false;  // Rmin is concave
2582        break;
2583      case kSPhi:
2584        if (fDPhi<=pi)     // Normal to Phi-
2585        {
2586          *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
2587          *validNorm=true;
2588        }
2589        else *validNorm=false;
2590        break ;
2591      case kEPhi:
2592        if (fDPhi<=pi)      // Normal to Phi+
2593        {
2594          *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
2595          *validNorm=true;
2596        }
2597        else *validNorm=false;
2598        break;
2599      case kSTheta:
2600        if( fSTheta == pi*0.5 )
2601        {
2602          *n=G4ThreeVector(0,0,1);
2603          *validNorm=true;
2604        }
2605        else if ( fSTheta > pi )
2606        {
2607          xi=p.x()+snxt*v.x();
2608          yi=p.y()+snxt*v.y();
2609          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2)) ;
2610          *n = G4ThreeVector(-xi/rhoSecTheta,   // N-
2611                             -yi/rhoSecTheta,
2612                             tanSTheta/std::sqrt(1+tanSTheta2)) ;
2613          *validNorm=true;
2614        }
2615        else *validNorm=false;  // Concave STheta cone
2616        break;
2617      case kETheta:
2618        if( ( fSTheta + fDTheta ) == pi*0.5 )
2619        {
2620          *n         = G4ThreeVector(0,0,-1);
2621          *validNorm = true ;
2622        }
2623        else if ( ( fSTheta + fDTheta ) < pi )
2624        {
2625          xi=p.x()+snxt*v.x();
2626          yi=p.y()+snxt*v.y();
2627          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2)) ;
2628          *n = G4ThreeVector( xi/rhoSecTheta,   // N+
2629                              yi/rhoSecTheta,
2630                              -tanSTheta/std::sqrt(1+tanSTheta2) ) ;
2631          *validNorm=true;
2632        }
2633        else *validNorm=false;   // Concave ETheta cone
2634        break;
2635      default:
2636        G4cout.precision(16);
2637        G4cout << G4endl;
2638        DumpInfo();
2639        G4cout << "Position:"  << G4endl << G4endl;
2640        G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
2641        G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
2642        G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
2643        G4cout << "Direction:" << G4endl << G4endl;
2644        G4cout << "v.x() = "   << v.x() << G4endl;
2645        G4cout << "v.y() = "   << v.y() << G4endl;
2646        G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
2647        G4cout << "Proposed distance :" << G4endl << G4endl;
2648        G4cout << "snxt = "    << snxt/mm << " mm" << G4endl << G4endl;
2649        G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2650                    "Notification", JustWarning,
2651                    "Undefined side for valid surface normal to solid.");
2652        break;
2653    }
2654  }
2655  if (snxt == kInfinity)
2656  {
2657    G4cout.precision(24);
2658    G4cout << G4endl;
2659    DumpInfo();
2660    G4cout << "Position:"  << G4endl << G4endl;
2661    G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
2662    G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
2663    G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
2664    G4cout << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm << " mm" 
2665           << G4endl << G4endl;
2666    G4cout << "Direction:" << G4endl << G4endl;
2667    G4cout << "v.x() = "   << v.x() << G4endl;
2668    G4cout << "v.y() = "   << v.y() << G4endl;
2669    G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
2670    G4cout << "Proposed distance :" << G4endl << G4endl;
2671    G4cout << "snxt = "    << snxt/mm << " mm" << G4endl << G4endl;
2672    G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2673                "Notification", JustWarning,
2674                "Logic error: snxt = kInfinity  ???");
2675  }
2676
2677  return snxt;
2678}
2679
2680/////////////////////////////////////////////////////////////////////////
2681//
2682// Calcluate distance (<=actual) to closest surface of shape from inside
2683
2684G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const
2685{
2686  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
2687  G4double rho2,rad,rho;
2688  G4double phiC,cosPhiC,sinPhiC,ePhi;
2689  G4double pTheta,dTheta1,dTheta2;
2690  rho2=p.x()*p.x()+p.y()*p.y();
2691  rad=std::sqrt(rho2+p.z()*p.z());
2692  rho=std::sqrt(rho2);
2693
2694#ifdef G4CSGDEBUG
2695  if( Inside(p) == kOutside )
2696  {
2697     G4cout.precision(16) ;
2698     G4cout << G4endl ;
2699     DumpInfo();
2700     G4cout << "Position:"  << G4endl << G4endl ;
2701     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
2702     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
2703     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
2704     G4Exception("G4Sphere::DistanceToOut(p)",
2705                 "Notification", JustWarning, "Point p is outside !?" );
2706  }
2707#endif
2708
2709  //
2710  // Distance to r shells
2711  //   
2712  if (fRmin)
2713  {
2714    safeRMin=rad-fRmin;
2715    safeRMax=fRmax-rad;
2716    if (safeRMin<safeRMax)
2717    {
2718      safe=safeRMin;
2719    }
2720    else
2721    {
2722      safe=safeRMax;
2723    }
2724  }
2725  else
2726  {
2727    safe=fRmax-rad;
2728  }
2729
2730  //
2731  // Distance to phi extent
2732  //
2733  if (fDPhi<twopi && rho)
2734  {
2735    phiC=fSPhi+fDPhi*0.5;
2736    cosPhiC=std::cos(phiC);
2737    sinPhiC=std::sin(phiC);
2738    if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
2739    {
2740      safePhi=-(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
2741    }
2742    else
2743    {
2744      ePhi=fSPhi+fDPhi;
2745      safePhi=(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
2746    }
2747    if (safePhi<safe) safe=safePhi;
2748  }
2749
2750  //
2751  // Distance to Theta extent
2752  //   
2753  if (rad)
2754  {
2755    pTheta=std::acos(p.z()/rad);
2756    if (pTheta<0) pTheta+=pi;
2757    dTheta1=pTheta-fSTheta;
2758    dTheta2=(fSTheta+fDTheta)-pTheta;
2759    if (dTheta1<dTheta2)
2760    {
2761      safeTheta=rad*std::sin(dTheta1);
2762      if (safe>safeTheta)
2763      {
2764        safe=safeTheta;
2765      }
2766    }
2767    else
2768    {
2769      safeTheta=rad*std::sin(dTheta2);
2770      if (safe>safeTheta)
2771      {
2772        safe=safeTheta;
2773      }
2774    }
2775  }
2776
2777  if (safe<0) safe=0;
2778    return safe;
2779}
2780
2781//////////////////////////////////////////////////////////////////////////
2782//
2783// Create a List containing the transformed vertices
2784// Ordering [0-3] -fDz cross section
2785//          [4-7] +fDz cross section such that [0] is below [4],
2786//                                             [1] below [5] etc.
2787// Note:
2788//  Caller has deletion resposibility
2789//  Potential improvement: For last slice, use actual ending angle
2790//                         to avoid rounding error problems.
2791
2792G4ThreeVectorList*
2793G4Sphere::CreateRotatedVertices( const G4AffineTransform& pTransform,
2794                                       G4int& noPolygonVertices ) const
2795{
2796  G4ThreeVectorList *vertices;
2797  G4ThreeVector vertex;
2798  G4double meshAnglePhi,meshRMax,crossAnglePhi,
2799           coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
2800  G4double meshTheta,crossTheta,startTheta;
2801  G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
2802  G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;
2803
2804  // Phi cross sections
2805   
2806  noPhiCrossSections=G4int (fDPhi/kMeshAngleDefault)+1;
2807   
2808  if (noPhiCrossSections<kMinMeshSections)
2809  {
2810    noPhiCrossSections=kMinMeshSections;
2811  }
2812  else if (noPhiCrossSections>kMaxMeshSections)
2813  {
2814    noPhiCrossSections=kMaxMeshSections;
2815  }
2816  meshAnglePhi=fDPhi/(noPhiCrossSections-1);
2817   
2818  // If complete in phi, set start angle such that mesh will be at fRMax
2819  // on the x axis. Will give better extent calculations when not rotated.
2820   
2821  if (fDPhi==pi*2.0 && fSPhi==0)
2822  {
2823    sAnglePhi = -meshAnglePhi*0.5;
2824  }
2825    else
2826  {
2827    sAnglePhi=fSPhi;
2828  }   
2829
2830  // Theta cross sections
2831   
2832  noThetaSections = G4int(fDTheta/kMeshAngleDefault)+1;
2833   
2834  if (noThetaSections<kMinMeshSections)
2835  {
2836    noThetaSections=kMinMeshSections;
2837  }
2838  else if (noThetaSections>kMaxMeshSections)
2839  {
2840    noThetaSections=kMaxMeshSections;
2841  }
2842  meshTheta=fDTheta/(noThetaSections-1);
2843   
2844  // If complete in Theta, set start angle such that mesh will be at fRMax
2845  // on the z axis. Will give better extent calculations when not rotated.
2846   
2847  if (fDTheta==pi && fSTheta==0)
2848  {
2849    startTheta = -meshTheta*0.5;
2850  }
2851  else
2852  {
2853    startTheta=fSTheta;
2854  }   
2855
2856  meshRMax = (meshAnglePhi >= meshTheta) ?
2857             fRmax/std::cos(meshAnglePhi*0.5) : fRmax/std::cos(meshTheta*0.5);
2858  G4double* cosCrossTheta = new G4double[noThetaSections];
2859  G4double* sinCrossTheta = new G4double[noThetaSections];   
2860  vertices=new G4ThreeVectorList();
2861  vertices->reserve(noPhiCrossSections*(noThetaSections*2));
2862  if (vertices && cosCrossTheta && sinCrossTheta)
2863  {
2864    for (crossSectionPhi=0;
2865         crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
2866    {
2867      crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
2868      coscrossAnglePhi=std::cos(crossAnglePhi);
2869      sincrossAnglePhi=std::sin(crossAnglePhi);
2870      for (crossSectionTheta=0;
2871           crossSectionTheta<noThetaSections;crossSectionTheta++)
2872      {
2873        // Compute coordinates of cross section at section crossSectionPhi
2874        //
2875        crossTheta=startTheta+crossSectionTheta*meshTheta;
2876        cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
2877        sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
2878
2879        rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
2880        rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
2881        rMinZ=fRmin*cosCrossTheta[crossSectionTheta];
2882       
2883        vertex=G4ThreeVector(rMinX,rMinY,rMinZ);
2884        vertices->push_back(pTransform.TransformPoint(vertex));
2885       
2886      }    // Theta forward
2887   
2888      for (crossSectionTheta=noThetaSections-1;
2889           crossSectionTheta>=0; crossSectionTheta--)
2890      {
2891        rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
2892        rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
2893        rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
2894       
2895        vertex=G4ThreeVector(rMaxX,rMaxY,rMaxZ);
2896        vertices->push_back(pTransform.TransformPoint(vertex));
2897
2898      }   // Theta back
2899    }   // Phi
2900    noPolygonVertices = noThetaSections*2 ;
2901  }
2902  else
2903  {
2904    DumpInfo();
2905    G4Exception("G4Sphere::CreateRotatedVertices()",
2906                "FatalError", FatalException,
2907                "Error in allocation of vertices. Out of memory !");
2908  }
2909
2910  delete[] cosCrossTheta;
2911  delete[] sinCrossTheta;
2912
2913  return vertices;
2914}
2915
2916//////////////////////////////////////////////////////////////////////////
2917//
2918// G4EntityType
2919
2920G4GeometryType G4Sphere::GetEntityType() const
2921{
2922  return G4String("G4Sphere");
2923}
2924
2925//////////////////////////////////////////////////////////////////////////
2926//
2927// Stream object contents to an output stream
2928
2929std::ostream& G4Sphere::StreamInfo( std::ostream& os ) const
2930{
2931  os << "-----------------------------------------------------------\n"
2932     << "    *** Dump for solid - " << GetName() << " ***\n"
2933     << "    ===================================================\n"
2934     << " Solid type: G4Sphere\n"
2935     << " Parameters: \n"
2936     << "    inner radius: " << fRmin/mm << " mm \n"
2937     << "    outer radius: " << fRmax/mm << " mm \n"
2938     << "    starting phi of segment  : " << fSPhi/degree << " degrees \n"
2939     << "    delta phi of segment     : " << fDPhi/degree << " degrees \n"
2940     << "    starting theta of segment: " << fSTheta/degree << " degrees \n"
2941     << "    delta theta of segment   : " << fDTheta/degree << " degrees \n"
2942     << "-----------------------------------------------------------\n";
2943
2944  return os;
2945}
2946
2947////////////////////////////////////////////////////////////////////////////////
2948//
2949// GetPointOnSurface
2950
2951G4ThreeVector G4Sphere::GetPointOnSurface() const
2952{
2953  G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
2954  G4double height1, height2, slant1, slant2, costheta, sintheta,theta,rRand;
2955
2956  height1 = (fRmax-fRmin)*std::cos(fSTheta);
2957  height2 = (fRmax-fRmin)*std::cos(fSTheta+fDTheta);
2958  slant1  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta))
2959                      + height1*height1);
2960  slant2  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta+fDTheta))
2961                      + height2*height2);
2962  rRand   = RandFlat::shoot(fRmin,fRmax);
2963 
2964  aOne = fRmax*fRmax*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
2965  aTwo = fRmin*fRmin*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
2966  aThr = fDPhi*((fRmax + fRmin)*std::sin(fSTheta))*slant1;
2967  aFou = fDPhi*((fRmax + fRmin)*std::sin(fSTheta+fDTheta))*slant2;
2968  aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
2969 
2970  phi = RandFlat::shoot(fSPhi, fSPhi + fDPhi); 
2971  cosphi = std::cos(phi); 
2972  sinphi = std::sin(phi);
2973  theta = RandFlat::shoot(fSTheta,fSTheta+fDTheta);
2974  costheta = std::cos(theta);
2975  sintheta = std::sqrt(1.-sqr(costheta));
2976
2977  if( ((fSPhi==0) && (fDPhi==2.*pi)) || (fDPhi==2.*pi) ) {aFiv = 0;}
2978  if(fSTheta == 0)  {aThr=0;}
2979  if(fDTheta + fSTheta == pi) {aFou = 0;}
2980  if(fSTheta == 0.5*pi) {aThr = pi*(fRmax*fRmax-fRmin*fRmin);}
2981  if(fSTheta + fDTheta == 0.5*pi) { aFou = pi*(fRmax*fRmax-fRmin*fRmin);}
2982
2983  chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
2984  if( (chose>=0.) && (chose<aOne) )
2985  {
2986    return G4ThreeVector(fRmax*sintheta*cosphi,
2987                         fRmax*sintheta*sinphi, fRmax*costheta);
2988  }
2989  else if( (chose>=aOne) && (chose<aOne+aTwo) )
2990  {
2991    return G4ThreeVector(fRmin*sintheta*cosphi,
2992                         fRmin*sintheta*sinphi, fRmin*costheta);
2993  }
2994  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
2995  {
2996    if (fSTheta != 0.5*pi)
2997    {
2998      zRand = RandFlat::shoot(fRmin*std::cos(fSTheta),fRmax*std::cos(fSTheta));
2999      return G4ThreeVector(std::tan(fSTheta)*zRand*cosphi,
3000                           std::tan(fSTheta)*zRand*sinphi,zRand);
3001    }
3002    else
3003    {
3004      return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
3005    }   
3006  }
3007  else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
3008  {
3009    if(fSTheta + fDTheta != 0.5*pi)
3010    {
3011      zRand = RandFlat::shoot(fRmin*std::cos(fSTheta+fDTheta),
3012                              fRmax*std::cos(fSTheta+fDTheta));
3013      return G4ThreeVector  (std::tan(fSTheta+fDTheta)*zRand*cosphi,
3014                             std::tan(fSTheta+fDTheta)*zRand*sinphi,zRand);
3015    }
3016    else
3017    {
3018      return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
3019    }
3020  }
3021  else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
3022  {
3023    return G4ThreeVector(rRand*sintheta*std::cos(fSPhi),
3024                         rRand*sintheta*std::sin(fSPhi),rRand*costheta);
3025  }
3026  else
3027  {
3028    return G4ThreeVector(rRand*sintheta*std::cos(fSPhi+fDPhi),
3029                         rRand*sintheta*std::sin(fSPhi+fDPhi),rRand*costheta);
3030  }
3031}
3032
3033/////////////////////////////////////////////////////////////////////////////
3034//
3035// Methods for visualisation
3036
3037G4VisExtent G4Sphere::GetExtent() const
3038{
3039  return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
3040}
3041
3042
3043void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
3044{
3045  scene.AddSolid (*this);
3046}
3047
3048G4Polyhedron* G4Sphere::CreatePolyhedron () const
3049{
3050  return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
3051}
3052
3053G4NURBS* G4Sphere::CreateNURBS () const
3054{
3055  return new G4NURBSbox (fRmax, fRmax, fRmax);       // Box for now!!!
3056}
Note: See TracBrowser for help on using the repository browser.