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