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