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