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