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