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