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

Last change on this file was 1337, checked in by garnier, 14 years ago

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