Geant4  10.02.p03
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 102294 2017-01-20 11:41:52Z gcosmo $
28 //
29 // class G4Sphere
30 //
31 // Implementation for G4Sphere class
32 //
33 // History:
34 //
35 // 05.04.12 M.Kelsey: GetPointOnSurface() throw flat in cos(theta), sqrt(r)
36 // 14.09.09 T.Nikitina: fix for phi section in DistanceToOut(p,v,..),as for G4Tubs,G4Cons
37 // 26.03.09 G.Cosmo : optimisations and uniform use of local radial tolerance
38 // 12.06.08 V.Grichine: fix for theta intersections in DistanceToOut(p,v,...)
39 // 22.07.05 O.Link : Added check for intersection with double cone
40 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
41 // 16.09.04 V.Grichine: bug fixed in SurfaceNormal(p), theta normals
42 // 16.07.04 V.Grichine: bug fixed in DistanceToOut(p,v), Rmin go outside
43 // 02.06.04 V.Grichine: bug fixed in DistanceToIn(p,v), on Rmax,Rmin go inside
44 // 30.10.03 J.Apostolakis: new algorithm in Inside for SPhi-sections
45 // 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi<SPhi, SPhi<0
46 // 19.06.02 V.Grichine: bug fixed in Inside(p), && -> && fDTheta - kAngTolerance
47 // 30.01.02 V.Grichine: bug fixed in Inside(p), && -> || at l.451
48 // 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
49 // 18.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
50 // 25.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), phi intersections
51 // 12.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), theta intersections
52 // 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
53 // 17.09.96 V.Grichine: final modifications to commit
54 // 28.03.94 P.Kent: old C++ code converted to tolerant geometry
55 // --------------------------------------------------------------------
56 
57 #include "G4Sphere.hh"
58 
59 #if !defined(G4GEOM_USE_USPHERE)
60 
61 #include "G4VoxelLimits.hh"
62 #include "G4AffineTransform.hh"
63 #include "G4GeometryTolerance.hh"
64 
65 #include "G4VPVParameterisation.hh"
66 
67 #include "Randomize.hh"
68 
69 #include "meshdefs.hh"
70 
71 #include "G4VGraphicsScene.hh"
72 #include "G4VisExtent.hh"
73 
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  // Special case rad2 = 0 comparing with direction
941  //
942  if ((rad2!=0.0) || (fRmin!=0.0))
943  {
944  // Keep going for computation of distance...
945  }
946  else // Positioned on the sphere's origin
947  {
948  G4double vTheta = std::atan2(std::sqrt(v.x()*v.x()+v.y()*v.y()),v.z()) ;
949  if ( (vTheta < tolSTheta) || (vTheta > tolETheta) )
950  {
951  return snxt ; // kInfinity
952  }
953  return snxt = 0.0 ;
954  }
955  }
956 
957  // Outer spherical shell intersection
958  // - Only if outside tolerant fRmax
959  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
960  // - No intersect -> no intersection with G4Sphere
961  //
962  // Shell eqn: x^2+y^2+z^2=RSPH^2
963  //
964  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
965  //
966  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
967  // => rad2 +2sd(pDotV3d) +sd^2 =R^2
968  //
969  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
970 
971  c = rad2 - fRmax*fRmax ;
972 
973  if (c > fRmaxTolerance*fRmax)
974  {
975  // If outside tolerant boundary of outer G4Sphere
976  // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]
977 
978  d2 = pDotV3d*pDotV3d - c ;
979 
980  if ( d2 >= 0 )
981  {
982  sd = -pDotV3d - std::sqrt(d2) ;
983 
984  if (sd >= 0 )
985  {
986  if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen on
987  { // 64 bits systems. Split long distances and recompute
988  G4double fTerm = sd-std::fmod(sd,dRmax);
989  sd = fTerm + DistanceToIn(p+fTerm*v,v);
990  }
991  xi = p.x() + sd*v.x() ;
992  yi = p.y() + sd*v.y() ;
993  rhoi = std::sqrt(xi*xi + yi*yi) ;
994 
995  if (!fFullPhiSphere && rhoi) // Check phi intersection
996  {
997  cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
998 
999  if (cosPsi >= cosHDPhiOT)
1000  {
1001  if (!fFullThetaSphere) // Check theta intersection
1002  {
1003  zi = p.z() + sd*v.z() ;
1004 
1005  // rhoi & zi can never both be 0
1006  // (=>intersect at origin =>fRmax=0)
1007  //
1008  iTheta = std::atan2(rhoi,zi) ;
1009  if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1010  {
1011  return snxt = sd ;
1012  }
1013  }
1014  else
1015  {
1016  return snxt=sd;
1017  }
1018  }
1019  }
1020  else
1021  {
1022  if (!fFullThetaSphere) // Check theta intersection
1023  {
1024  zi = p.z() + sd*v.z() ;
1025 
1026  // rhoi & zi can never both be 0
1027  // (=>intersect at origin => fRmax=0 !)
1028  //
1029  iTheta = std::atan2(rhoi,zi) ;
1030  if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1031  {
1032  return snxt=sd;
1033  }
1034  }
1035  else
1036  {
1037  return snxt = sd;
1038  }
1039  }
1040  }
1041  }
1042  else // No intersection with G4Sphere
1043  {
1044  return snxt=kInfinity;
1045  }
1046  }
1047  else
1048  {
1049  // Inside outer radius
1050  // check not inside, and heading through G4Sphere (-> 0 to in)
1051 
1052  d2 = pDotV3d*pDotV3d - c ;
1053 
1054  if ( (rad2 > tolIRMax2)
1055  && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
1056  {
1057  if (!fFullPhiSphere)
1058  {
1059  // Use inner phi tolerant boundary -> if on tolerant
1060  // phi boundaries, phi intersect code handles leaving/entering checks
1061 
1062  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1063 
1064  if (cosPsi>=cosHDPhiIT)
1065  {
1066  // inside radii, delta r -ve, inside phi
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  else
1083  {
1084  if ( !fFullThetaSphere )
1085  {
1086  if ( (pTheta >= tolSTheta + kAngTolerance)
1087  && (pTheta <= tolETheta - kAngTolerance) )
1088  {
1089  return snxt=0;
1090  }
1091  }
1092  else // strictly inside Theta in both cases
1093  {
1094  return snxt=0;
1095  }
1096  }
1097  }
1098  }
1099 
1100  // Inner spherical shell intersection
1101  // - Always farthest root, because would have passed through outer
1102  // surface first.
1103  // - Tolerant check if travelling through solid
1104 
1105  if (fRmin)
1106  {
1107  c = rad2 - fRmin*fRmin ;
1108  d2 = pDotV3d*pDotV3d - c ;
1109 
1110  // Within tolerance inner radius of inner G4Sphere
1111  // Check for immediate entry/already inside and travelling outwards
1112 
1113  if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
1114  && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
1115  {
1116  if ( !fFullPhiSphere )
1117  {
1118  // Use inner phi tolerant boundary -> if on tolerant
1119  // phi boundaries, phi intersect code handles leaving/entering checks
1120 
1121  cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
1122  if (cosPsi >= cosHDPhiIT)
1123  {
1124  // inside radii, delta r -ve, inside phi
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
1141  {
1142  if ( !fFullThetaSphere )
1143  {
1144  if ( (pTheta >= tolSTheta + kAngTolerance)
1145  && (pTheta <= tolETheta - kAngTolerance) )
1146  {
1147  return snxt = 0 ;
1148  }
1149  }
1150  else
1151  {
1152  return snxt=0;
1153  }
1154  }
1155  }
1156  else // Not special tolerant case
1157  {
1158  if (d2 >= 0)
1159  {
1160  sd = -pDotV3d + std::sqrt(d2) ;
1161  if ( sd >= halfRminTolerance ) // It was >= 0 ??
1162  {
1163  xi = p.x() + sd*v.x() ;
1164  yi = p.y() + sd*v.y() ;
1165  rhoi = std::sqrt(xi*xi+yi*yi) ;
1166 
1167  if ( !fFullPhiSphere && rhoi ) // Check phi intersection
1168  {
1169  cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
1170 
1171  if (cosPsi >= cosHDPhiOT)
1172  {
1173  if ( !fFullThetaSphere ) // Check theta intersection
1174  {
1175  zi = p.z() + sd*v.z() ;
1176 
1177  // rhoi & zi can never both be 0
1178  // (=>intersect at origin =>fRmax=0)
1179  //
1180  iTheta = std::atan2(rhoi,zi) ;
1181  if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
1182  {
1183  snxt = sd;
1184  }
1185  }
1186  else
1187  {
1188  snxt=sd;
1189  }
1190  }
1191  }
1192  else
1193  {
1194  if ( !fFullThetaSphere ) // Check theta intersection
1195  {
1196  zi = p.z() + sd*v.z() ;
1197 
1198  // rhoi & zi can never both be 0
1199  // (=>intersect at origin => fRmax=0 !)
1200  //
1201  iTheta = std::atan2(rhoi,zi) ;
1202  if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1203  {
1204  snxt = sd;
1205  }
1206  }
1207  else
1208  {
1209  snxt = sd;
1210  }
1211  }
1212  }
1213  }
1214  }
1215  }
1216 
1217  // Phi segment intersection
1218  //
1219  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1220  //
1221  // o NOTE: Large duplication of code between sphi & ephi checks
1222  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1223  // intersection check <=0 -> >=0
1224  // -> Should use some form of loop Construct
1225  //
1226  if ( !fFullPhiSphere )
1227  {
1228  // First phi surface ('S'tarting phi)
1229  // Comp = Component in outwards normal dirn
1230  //
1231  Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1232 
1233  if ( Comp < 0 )
1234  {
1235  Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
1236 
1237  if (Dist < halfCarTolerance)
1238  {
1239  sd = Dist/Comp ;
1240 
1241  if (sd < snxt)
1242  {
1243  if ( sd > 0 )
1244  {
1245  xi = p.x() + sd*v.x() ;
1246  yi = p.y() + sd*v.y() ;
1247  zi = p.z() + sd*v.z() ;
1248  rhoi2 = xi*xi + yi*yi ;
1249  radi2 = rhoi2 + zi*zi ;
1250  }
1251  else
1252  {
1253  sd = 0 ;
1254  xi = p.x() ;
1255  yi = p.y() ;
1256  zi = p.z() ;
1257  rhoi2 = rho2 ;
1258  radi2 = rad2 ;
1259  }
1260  if ( (radi2 <= tolORMax2)
1261  && (radi2 >= tolORMin2)
1262  && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
1263  {
1264  // Check theta intersection
1265  // rhoi & zi can never both be 0
1266  // (=>intersect at origin =>fRmax=0)
1267  //
1268  if ( !fFullThetaSphere )
1269  {
1270  iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1271  if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1272  {
1273  // r and theta intersections good
1274  // - check intersecting with correct half-plane
1275 
1276  if ((yi*cosCPhi-xi*sinCPhi) <= 0)
1277  {
1278  snxt = sd;
1279  }
1280  }
1281  }
1282  else
1283  {
1284  snxt = sd;
1285  }
1286  }
1287  }
1288  }
1289  }
1290 
1291  // Second phi surface ('E'nding phi)
1292  // Component in outwards normal dirn
1293 
1294  Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
1295 
1296  if (Comp < 0)
1297  {
1298  Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
1299  if ( Dist < halfCarTolerance )
1300  {
1301  sd = Dist/Comp ;
1302 
1303  if ( sd < snxt )
1304  {
1305  if (sd > 0)
1306  {
1307  xi = p.x() + sd*v.x() ;
1308  yi = p.y() + sd*v.y() ;
1309  zi = p.z() + sd*v.z() ;
1310  rhoi2 = xi*xi + yi*yi ;
1311  radi2 = rhoi2 + zi*zi ;
1312  }
1313  else
1314  {
1315  sd = 0 ;
1316  xi = p.x() ;
1317  yi = p.y() ;
1318  zi = p.z() ;
1319  rhoi2 = rho2 ;
1320  radi2 = rad2 ;
1321  }
1322  if ( (radi2 <= tolORMax2)
1323  && (radi2 >= tolORMin2)
1324  && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
1325  {
1326  // Check theta intersection
1327  // rhoi & zi can never both be 0
1328  // (=>intersect at origin =>fRmax=0)
1329  //
1330  if ( !fFullThetaSphere )
1331  {
1332  iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
1333  if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
1334  {
1335  // r and theta intersections good
1336  // - check intersecting with correct half-plane
1337 
1338  if ((yi*cosCPhi-xi*sinCPhi) >= 0)
1339  {
1340  snxt = sd;
1341  }
1342  }
1343  }
1344  else
1345  {
1346  snxt = sd;
1347  }
1348  }
1349  }
1350  }
1351  }
1352  }
1353 
1354  // Theta segment intersection
1355 
1356  if ( !fFullThetaSphere )
1357  {
1358 
1359  // Intersection with theta surfaces
1360  // Known failure cases:
1361  // o Inside tolerance of stheta surface, skim
1362  // ~parallel to cone and Hit & enter etheta surface [& visa versa]
1363  //
1364  // To solve: Check 2nd root of etheta surface in addition to stheta
1365  //
1366  // o start/end theta is exactly pi/2
1367  // Intersections with cones
1368  //
1369  // Cone equation: x^2+y^2=z^2tan^2(t)
1370  //
1371  // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
1372  //
1373  // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
1374  // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
1375  //
1376  // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
1377 
1378  if (fSTheta)
1379  {
1380  dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
1381  }
1382  else
1383  {
1384  dist2STheta = kInfinity ;
1385  }
1386  if ( eTheta < pi )
1387  {
1388  dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
1389  }
1390  else
1391  {
1392  dist2ETheta=kInfinity;
1393  }
1394  if ( pTheta < tolSTheta )
1395  {
1396  // Inside (theta<stheta-tol) stheta cone
1397  // First root of stheta cone, second if first root -ve
1398 
1399  t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1400  t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1401  if (t1)
1402  {
1403  b = t2/t1 ;
1404  c = dist2STheta/t1 ;
1405  d2 = b*b - c ;
1406 
1407  if ( d2 >= 0 )
1408  {
1409  d = std::sqrt(d2) ;
1410  sd = -b - d ; // First root
1411  zi = p.z() + sd*v.z();
1412 
1413  if ( (sd < 0) || (zi*(fSTheta - halfpi) > 0) )
1414  {
1415  sd = -b+d; // Second root
1416  }
1417  if ((sd >= 0) && (sd < snxt))
1418  {
1419  xi = p.x() + sd*v.x();
1420  yi = p.y() + sd*v.y();
1421  zi = p.z() + sd*v.z();
1422  rhoi2 = xi*xi + yi*yi;
1423  radi2 = rhoi2 + zi*zi;
1424  if ( (radi2 <= tolORMax2)
1425  && (radi2 >= tolORMin2)
1426  && (zi*(fSTheta - halfpi) <= 0) )
1427  {
1428  if ( !fFullPhiSphere && rhoi2 ) // Check phi intersection
1429  {
1430  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1431  if (cosPsi >= cosHDPhiOT)
1432  {
1433  snxt = sd;
1434  }
1435  }
1436  else
1437  {
1438  snxt = sd;
1439  }
1440  }
1441  }
1442  }
1443  }
1444 
1445  // Possible intersection with ETheta cone.
1446  // Second >= 0 root should be considered
1447 
1448  if ( eTheta < pi )
1449  {
1450  t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1451  t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1452  if (t1)
1453  {
1454  b = t2/t1 ;
1455  c = dist2ETheta/t1 ;
1456  d2 = b*b - c ;
1457 
1458  if (d2 >= 0)
1459  {
1460  d = std::sqrt(d2) ;
1461  sd = -b + d ; // Second root
1462 
1463  if ( (sd >= 0) && (sd < snxt) )
1464  {
1465  xi = p.x() + sd*v.x() ;
1466  yi = p.y() + sd*v.y() ;
1467  zi = p.z() + sd*v.z() ;
1468  rhoi2 = xi*xi + yi*yi ;
1469  radi2 = rhoi2 + zi*zi ;
1470 
1471  if ( (radi2 <= tolORMax2)
1472  && (radi2 >= tolORMin2)
1473  && (zi*(eTheta - halfpi) <= 0) )
1474  {
1475  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1476  {
1477  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1478  if (cosPsi >= cosHDPhiOT)
1479  {
1480  snxt = sd;
1481  }
1482  }
1483  else
1484  {
1485  snxt = sd;
1486  }
1487  }
1488  }
1489  }
1490  }
1491  }
1492  }
1493  else if ( pTheta > tolETheta )
1494  {
1495  // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
1496  // Inside (theta > etheta+tol) e-theta cone
1497  // First root of etheta cone, second if first root 'imaginary'
1498 
1499  t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1500  t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1501  if (t1)
1502  {
1503  b = t2/t1 ;
1504  c = dist2ETheta/t1 ;
1505  d2 = b*b - c ;
1506 
1507  if (d2 >= 0)
1508  {
1509  d = std::sqrt(d2) ;
1510  sd = -b - d ; // First root
1511  zi = p.z() + sd*v.z();
1512 
1513  if ( (sd < 0) || (zi*(eTheta - halfpi) > 0) )
1514  {
1515  sd = -b + d ; // second root
1516  }
1517  if ( (sd >= 0) && (sd < snxt) )
1518  {
1519  xi = p.x() + sd*v.x() ;
1520  yi = p.y() + sd*v.y() ;
1521  zi = p.z() + sd*v.z() ;
1522  rhoi2 = xi*xi + yi*yi ;
1523  radi2 = rhoi2 + zi*zi ;
1524 
1525  if ( (radi2 <= tolORMax2)
1526  && (radi2 >= tolORMin2)
1527  && (zi*(eTheta - halfpi) <= 0) )
1528  {
1529  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1530  {
1531  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1532  if (cosPsi >= cosHDPhiOT)
1533  {
1534  snxt = sd;
1535  }
1536  }
1537  else
1538  {
1539  snxt = sd;
1540  }
1541  }
1542  }
1543  }
1544  }
1545 
1546  // Possible intersection with STheta cone.
1547  // Second >= 0 root should be considered
1548 
1549  if ( fSTheta )
1550  {
1551  t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1552  t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1553  if (t1)
1554  {
1555  b = t2/t1 ;
1556  c = dist2STheta/t1 ;
1557  d2 = b*b - c ;
1558 
1559  if (d2 >= 0)
1560  {
1561  d = std::sqrt(d2) ;
1562  sd = -b + d ; // Second root
1563 
1564  if ( (sd >= 0) && (sd < snxt) )
1565  {
1566  xi = p.x() + sd*v.x() ;
1567  yi = p.y() + sd*v.y() ;
1568  zi = p.z() + sd*v.z() ;
1569  rhoi2 = xi*xi + yi*yi ;
1570  radi2 = rhoi2 + zi*zi ;
1571 
1572  if ( (radi2 <= tolORMax2)
1573  && (radi2 >= tolORMin2)
1574  && (zi*(fSTheta - halfpi) <= 0) )
1575  {
1576  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1577  {
1578  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1579  if (cosPsi >= cosHDPhiOT)
1580  {
1581  snxt = sd;
1582  }
1583  }
1584  else
1585  {
1586  snxt = sd;
1587  }
1588  }
1589  }
1590  }
1591  }
1592  }
1593  }
1594  else if ( (pTheta < tolSTheta + kAngTolerance)
1595  && (fSTheta > halfAngTolerance) )
1596  {
1597  // In tolerance of stheta
1598  // If entering through solid [r,phi] => 0 to in
1599  // else try 2nd root
1600 
1601  t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1602  if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
1603  || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
1604  || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
1605  {
1606  if (!fFullPhiSphere && rho2) // Check phi intersection
1607  {
1608  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1609  if (cosPsi >= cosHDPhiIT)
1610  {
1611  return 0 ;
1612  }
1613  }
1614  else
1615  {
1616  return 0 ;
1617  }
1618  }
1619 
1620  // Not entering immediately/travelling through
1621 
1622  t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1623  if (t1)
1624  {
1625  b = t2/t1 ;
1626  c = dist2STheta/t1 ;
1627  d2 = b*b - c ;
1628 
1629  if (d2 >= 0)
1630  {
1631  d = std::sqrt(d2) ;
1632  sd = -b + d ;
1633  if ( (sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < halfpi) )
1634  { // ^^^^^^^^^^^^^^^^^^^^^ shouldn't it be >=0 instead ?
1635  xi = p.x() + sd*v.x() ;
1636  yi = p.y() + sd*v.y() ;
1637  zi = p.z() + sd*v.z() ;
1638  rhoi2 = xi*xi + yi*yi ;
1639  radi2 = rhoi2 + zi*zi ;
1640 
1641  if ( (radi2 <= tolORMax2)
1642  && (radi2 >= tolORMin2)
1643  && (zi*(fSTheta - halfpi) <= 0) )
1644  {
1645  if ( !fFullPhiSphere && rhoi2 ) // Check phi intersection
1646  {
1647  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1648  if ( cosPsi >= cosHDPhiOT )
1649  {
1650  snxt = sd;
1651  }
1652  }
1653  else
1654  {
1655  snxt = sd;
1656  }
1657  }
1658  }
1659  }
1660  }
1661  }
1662  else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
1663  {
1664 
1665  // In tolerance of etheta
1666  // If entering through solid [r,phi] => 0 to in
1667  // else try 2nd root
1668 
1669  t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1670 
1671  if ( ((t2<0) && (eTheta < halfpi)
1672  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1673  || ((t2>=0) && (eTheta > halfpi)
1674  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
1675  || ((v.z()>0) && (eTheta == halfpi)
1676  && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)) )
1677  {
1678  if (!fFullPhiSphere && rho2) // Check phi intersection
1679  {
1680  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
1681  if (cosPsi >= cosHDPhiIT)
1682  {
1683  return 0 ;
1684  }
1685  }
1686  else
1687  {
1688  return 0 ;
1689  }
1690  }
1691 
1692  // Not entering immediately/travelling through
1693 
1694  t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1695  if (t1)
1696  {
1697  b = t2/t1 ;
1698  c = dist2ETheta/t1 ;
1699  d2 = b*b - c ;
1700 
1701  if (d2 >= 0)
1702  {
1703  d = std::sqrt(d2) ;
1704  sd = -b + d ;
1705 
1706  if ( (sd >= halfCarTolerance)
1707  && (sd < snxt) && (eTheta > halfpi) )
1708  {
1709  xi = p.x() + sd*v.x() ;
1710  yi = p.y() + sd*v.y() ;
1711  zi = p.z() + sd*v.z() ;
1712  rhoi2 = xi*xi + yi*yi ;
1713  radi2 = rhoi2 + zi*zi ;
1714 
1715  if ( (radi2 <= tolORMax2)
1716  && (radi2 >= tolORMin2)
1717  && (zi*(eTheta - halfpi) <= 0) )
1718  {
1719  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1720  {
1721  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1722  if (cosPsi >= cosHDPhiOT)
1723  {
1724  snxt = sd;
1725  }
1726  }
1727  else
1728  {
1729  snxt = sd;
1730  }
1731  }
1732  }
1733  }
1734  }
1735  }
1736  else
1737  {
1738  // stheta+tol<theta<etheta-tol
1739  // For BOTH stheta & etheta check 2nd root for validity [r,phi]
1740 
1741  t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
1742  t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
1743  if (t1)
1744  {
1745  b = t2/t1;
1746  c = dist2STheta/t1 ;
1747  d2 = b*b - c ;
1748 
1749  if (d2 >= 0)
1750  {
1751  d = std::sqrt(d2) ;
1752  sd = -b + d ; // second root
1753 
1754  if ((sd >= 0) && (sd < snxt))
1755  {
1756  xi = p.x() + sd*v.x() ;
1757  yi = p.y() + sd*v.y() ;
1758  zi = p.z() + sd*v.z() ;
1759  rhoi2 = xi*xi + yi*yi ;
1760  radi2 = rhoi2 + zi*zi ;
1761 
1762  if ( (radi2 <= tolORMax2)
1763  && (radi2 >= tolORMin2)
1764  && (zi*(fSTheta - halfpi) <= 0) )
1765  {
1766  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1767  {
1768  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1769  if (cosPsi >= cosHDPhiOT)
1770  {
1771  snxt = sd;
1772  }
1773  }
1774  else
1775  {
1776  snxt = sd;
1777  }
1778  }
1779  }
1780  }
1781  }
1782  t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
1783  t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
1784  if (t1)
1785  {
1786  b = t2/t1 ;
1787  c = dist2ETheta/t1 ;
1788  d2 = b*b - c ;
1789 
1790  if (d2 >= 0)
1791  {
1792  d = std::sqrt(d2) ;
1793  sd = -b + d; // second root
1794 
1795  if ((sd >= 0) && (sd < snxt))
1796  {
1797  xi = p.x() + sd*v.x() ;
1798  yi = p.y() + sd*v.y() ;
1799  zi = p.z() + sd*v.z() ;
1800  rhoi2 = xi*xi + yi*yi ;
1801  radi2 = rhoi2 + zi*zi ;
1802 
1803  if ( (radi2 <= tolORMax2)
1804  && (radi2 >= tolORMin2)
1805  && (zi*(eTheta - halfpi) <= 0) )
1806  {
1807  if (!fFullPhiSphere && rhoi2) // Check phi intersection
1808  {
1809  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
1810  if ( cosPsi >= cosHDPhiOT )
1811  {
1812  snxt = sd;
1813  }
1814  }
1815  else
1816  {
1817  snxt = sd;
1818  }
1819  }
1820  }
1821  }
1822  }
1823  }
1824  }
1825  return snxt;
1826 }
1827 
1829 //
1830 // Calculate distance (<= actual) to closest surface of shape from outside
1831 // - Calculate distance to radial planes
1832 // - Only to phi planes if outside phi extent
1833 // - Only to theta planes if outside theta extent
1834 // - Return 0 if point inside
1835 
1836 G4double G4Sphere::DistanceToIn( const G4ThreeVector& p ) const
1837 {
1838  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
1839  G4double rho2,rds,rho;
1840  G4double cosPsi;
1841  G4double pTheta,dTheta1,dTheta2;
1842  rho2=p.x()*p.x()+p.y()*p.y();
1843  rds=std::sqrt(rho2+p.z()*p.z());
1844  rho=std::sqrt(rho2);
1845 
1846  //
1847  // Distance to r shells
1848  //
1849  if (fRmin)
1850  {
1851  safeRMin=fRmin-rds;
1852  safeRMax=rds-fRmax;
1853  if (safeRMin>safeRMax)
1854  {
1855  safe=safeRMin;
1856  }
1857  else
1858  {
1859  safe=safeRMax;
1860  }
1861  }
1862  else
1863  {
1864  safe=rds-fRmax;
1865  }
1866 
1867  //
1868  // Distance to phi extent
1869  //
1870  if (!fFullPhiSphere && rho)
1871  {
1872  // Psi=angle from central phi to point
1873  //
1874  cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
1875  if (cosPsi<std::cos(hDPhi))
1876  {
1877  // Point lies outside phi range
1878  //
1879  if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
1880  {
1881  safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
1882  }
1883  else
1884  {
1885  safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1886  }
1887  if (safePhi>safe) { safe=safePhi; }
1888  }
1889  }
1890  //
1891  // Distance to Theta extent
1892  //
1893  if ((rds!=0.0) && (!fFullThetaSphere))
1894  {
1895  pTheta=std::acos(p.z()/rds);
1896  if (pTheta<0) { pTheta+=pi; }
1897  dTheta1=fSTheta-pTheta;
1898  dTheta2=pTheta-eTheta;
1899  if (dTheta1>dTheta2)
1900  {
1901  if (dTheta1>=0) // WHY ???????????
1902  {
1903  safeTheta=rds*std::sin(dTheta1);
1904  if (safe<=safeTheta)
1905  {
1906  safe=safeTheta;
1907  }
1908  }
1909  }
1910  else
1911  {
1912  if (dTheta2>=0)
1913  {
1914  safeTheta=rds*std::sin(dTheta2);
1915  if (safe<=safeTheta)
1916  {
1917  safe=safeTheta;
1918  }
1919  }
1920  }
1921  }
1922 
1923  if (safe<0) { safe=0; }
1924  return safe;
1925 }
1926 
1928 //
1929 // Calculate distance to surface of shape from 'inside', allowing for tolerance
1930 // - Only Calc rmax intersection if no valid rmin intersection
1931 
1932 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p,
1933  const G4ThreeVector& v,
1934  const G4bool calcNorm,
1935  G4bool *validNorm,
1936  G4ThreeVector *n ) const
1937 {
1938  G4double snxt = kInfinity; // snxt is default return value
1939  G4double sphi= kInfinity,stheta= kInfinity;
1940  ESide side=kNull,sidephi=kNull,sidetheta=kNull;
1941 
1942  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
1943  const G4double halfRminTolerance = fRminTolerance*0.5;
1944  const G4double Rmax_plus = fRmax + halfRmaxTolerance;
1945  const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
1946  G4double t1,t2;
1947  G4double b,c,d;
1948 
1949  // Variables for phi intersection:
1950 
1951  G4double pDistS,compS,pDistE,compE,sphi2,vphi;
1952 
1953  G4double rho2,rad2,pDotV2d,pDotV3d;
1954 
1955  G4double xi,yi,zi; // Intersection point
1956 
1957  // Theta precals
1958  //
1959  G4double rhoSecTheta;
1960  G4double dist2STheta, dist2ETheta, distTheta;
1961  G4double d2,sd;
1962 
1963  // General Precalcs
1964  //
1965  rho2 = p.x()*p.x()+p.y()*p.y();
1966  rad2 = rho2+p.z()*p.z();
1967 
1968  pDotV2d = p.x()*v.x()+p.y()*v.y();
1969  pDotV3d = pDotV2d+p.z()*v.z();
1970 
1971  // Radial Intersections from G4Sphere::DistanceToIn
1972  //
1973  // Outer spherical shell intersection
1974  // - Only if outside tolerant fRmax
1975  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
1976  // - No intersect -> no intersection with G4Sphere
1977  //
1978  // Shell eqn: x^2+y^2+z^2=RSPH^2
1979  //
1980  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
1981  //
1982  // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2
1983  // => rad2 +2sd(pDotV3d) +sd^2 =R^2
1984  //
1985  // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
1986 
1987  if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
1988  {
1989  c = rad2 - fRmax*fRmax;
1990 
1991  if (c < fRmaxTolerance*fRmax)
1992  {
1993  // Within tolerant Outer radius
1994  //
1995  // The test is
1996  // rad - fRmax < 0.5*kRadTolerance
1997  // => rad < fRmax + 0.5*kRadTol
1998  // => rad2 < (fRmax + 0.5*kRadTol)^2
1999  // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
2000  // => rad2 - fRmax^2 <~ fRmax*kRadTol
2001 
2002  d2 = pDotV3d*pDotV3d - c;
2003 
2004  if( (c >- fRmaxTolerance*fRmax) // on tolerant surface
2005  && ((pDotV3d >=0) || (d2 < 0)) ) // leaving outside from Rmax
2006  // not re-entering
2007  {
2008  if(calcNorm)
2009  {
2010  *validNorm = true ;
2011  *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
2012  }
2013  return snxt = 0;
2014  }
2015  else
2016  {
2017  snxt = -pDotV3d+std::sqrt(d2); // second root since inside Rmax
2018  side = kRMax ;
2019  }
2020  }
2021 
2022  // Inner spherical shell intersection:
2023  // Always first >=0 root, because would have passed
2024  // from outside of Rmin surface .
2025 
2026  if (fRmin)
2027  {
2028  c = rad2 - fRmin*fRmin;
2029  d2 = pDotV3d*pDotV3d - c;
2030 
2031  if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
2032  {
2033  if ( (c < fRminTolerance*fRmin) // leaving from Rmin
2034  && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
2035  {
2036  if(calcNorm) { *validNorm = false; } // Rmin surface is concave
2037  return snxt = 0 ;
2038  }
2039  else
2040  {
2041  if ( d2 >= 0. )
2042  {
2043  sd = -pDotV3d-std::sqrt(d2);
2044 
2045  if ( sd >= 0. ) // Always intersect Rmin first
2046  {
2047  snxt = sd ;
2048  side = kRMin ;
2049  }
2050  }
2051  }
2052  }
2053  }
2054  }
2055 
2056  // Theta segment intersection
2057 
2058  if ( !fFullThetaSphere )
2059  {
2060  // Intersection with theta surfaces
2061  //
2062  // Known failure cases:
2063  // o Inside tolerance of stheta surface, skim
2064  // ~parallel to cone and Hit & enter etheta surface [& visa versa]
2065  //
2066  // To solve: Check 2nd root of etheta surface in addition to stheta
2067  //
2068  // o start/end theta is exactly pi/2
2069  //
2070  // Intersections with cones
2071  //
2072  // Cone equation: x^2+y^2=z^2tan^2(t)
2073  //
2074  // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
2075  //
2076  // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t))
2077  // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0
2078  //
2079  // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
2080  //
2081 
2082  if(fSTheta) // intersection with first cons
2083  {
2084  if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
2085  {
2086  if( v.z() > 0. )
2087  {
2088  if ( std::fabs( p.z() ) <= halfRmaxTolerance )
2089  {
2090  if(calcNorm)
2091  {
2092  *validNorm = true;
2093  *n = G4ThreeVector(0.,0.,1.);
2094  }
2095  return snxt = 0 ;
2096  }
2097  stheta = -p.z()/v.z();
2098  sidetheta = kSTheta;
2099  }
2100  }
2101  else // kons is not plane
2102  {
2103  t1 = 1-v.z()*v.z()*(1+tanSTheta2);
2104  t2 = pDotV2d-p.z()*v.z()*tanSTheta2; // ~vDotN if p on cons
2105  dist2STheta = rho2-p.z()*p.z()*tanSTheta2; // t3
2106 
2107  distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
2108 
2109  if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
2110  { // v parallel to kons
2111  if( v.z() > 0. )
2112  {
2113  if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
2114  {
2115  if( (fSTheta < halfpi) && (p.z() > 0.) )
2116  {
2117  if( calcNorm ) { *validNorm = false; }
2118  return snxt = 0.;
2119  }
2120  else if( (fSTheta > halfpi) && (p.z() <= 0) )
2121  {
2122  if( calcNorm )
2123  {
2124  *validNorm = true;
2125  if (rho2)
2126  {
2127  rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
2128 
2129  *n = G4ThreeVector( p.x()/rhoSecTheta,
2130  p.y()/rhoSecTheta,
2131  std::sin(fSTheta) );
2132  }
2133  else *n = G4ThreeVector(0.,0.,1.);
2134  }
2135  return snxt = 0.;
2136  }
2137  }
2138  stheta = -0.5*dist2STheta/t2;
2139  sidetheta = kSTheta;
2140  }
2141  } // 2nd order equation, 1st root of fSTheta cone,
2142  else // 2nd if 1st root -ve
2143  {
2144  if( std::fabs(distTheta) < halfRmaxTolerance )
2145  {
2146  if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
2147  {
2148  if( calcNorm )
2149  {
2150  *validNorm = true;
2151  if (rho2)
2152  {
2153  rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
2154 
2155  *n = G4ThreeVector( p.x()/rhoSecTheta,
2156  p.y()/rhoSecTheta,
2157  std::sin(fSTheta) );
2158  }
2159  else { *n = G4ThreeVector(0.,0.,1.); }
2160  }
2161  return snxt = 0.;
2162  }
2163  else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
2164  {
2165  if( calcNorm ) { *validNorm = false; }
2166  return snxt = 0.;
2167  }
2168  }
2169  b = t2/t1;
2170  c = dist2STheta/t1;
2171  d2 = b*b - c ;
2172 
2173  if ( d2 >= 0. )
2174  {
2175  d = std::sqrt(d2);
2176 
2177  if( fSTheta > halfpi )
2178  {
2179  sd = -b - d; // First root
2180 
2181  if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
2182  || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
2183  {
2184  sd = -b + d ; // 2nd root
2185  }
2186  if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )
2187  {
2188  stheta = sd;
2189  sidetheta = kSTheta;
2190  }
2191  }
2192  else // sTheta < pi/2, concave surface, no normal
2193  {
2194  sd = -b - d; // First root
2195 
2196  if ( ( (std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.) )
2197  || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() < 0.) ) )
2198  {
2199  sd = -b + d ; // 2nd root
2200  }
2201  if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() >= 0.) )
2202  {
2203  stheta = sd;
2204  sidetheta = kSTheta;
2205  }
2206  }
2207  }
2208  }
2209  }
2210  }
2211  if (eTheta < pi) // intersection with second cons
2212  {
2213  if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
2214  {
2215  if( v.z() < 0. )
2216  {
2217  if ( std::fabs( p.z() ) <= halfRmaxTolerance )
2218  {
2219  if(calcNorm)
2220  {
2221  *validNorm = true;
2222  *n = G4ThreeVector(0.,0.,-1.);
2223  }
2224  return snxt = 0 ;
2225  }
2226  sd = -p.z()/v.z();
2227 
2228  if( sd < stheta )
2229  {
2230  stheta = sd;
2231  sidetheta = kETheta;
2232  }
2233  }
2234  }
2235  else // kons is not plane
2236  {
2237  t1 = 1-v.z()*v.z()*(1+tanETheta2);
2238  t2 = pDotV2d-p.z()*v.z()*tanETheta2; // ~vDotN if p on cons
2239  dist2ETheta = rho2-p.z()*p.z()*tanETheta2; // t3
2240 
2241  distTheta = std::sqrt(rho2)-p.z()*tanETheta;
2242 
2243  if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
2244  { // v parallel to kons
2245  if( v.z() < 0. )
2246  {
2247  if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
2248  {
2249  if( (eTheta > halfpi) && (p.z() < 0.) )
2250  {
2251  if( calcNorm ) { *validNorm = false; }
2252  return snxt = 0.;
2253  }
2254  else if ( (eTheta < halfpi) && (p.z() >= 0) )
2255  {
2256  if( calcNorm )
2257  {
2258  *validNorm = true;
2259  if (rho2)
2260  {
2261  rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2262  *n = G4ThreeVector( p.x()/rhoSecTheta,
2263  p.y()/rhoSecTheta,
2264  -sinETheta );
2265  }
2266  else { *n = G4ThreeVector(0.,0.,-1.); }
2267  }
2268  return snxt = 0.;
2269  }
2270  }
2271  sd = -0.5*dist2ETheta/t2;
2272 
2273  if( sd < stheta )
2274  {
2275  stheta = sd;
2276  sidetheta = kETheta;
2277  }
2278  }
2279  } // 2nd order equation, 1st root of fSTheta cone
2280  else // 2nd if 1st root -ve
2281  {
2282  if ( std::fabs(distTheta) < halfRmaxTolerance )
2283  {
2284  if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
2285  {
2286  if( calcNorm )
2287  {
2288  *validNorm = true;
2289  if (rho2)
2290  {
2291  rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2292  *n = G4ThreeVector( p.x()/rhoSecTheta,
2293  p.y()/rhoSecTheta,
2294  -sinETheta );
2295  }
2296  else *n = G4ThreeVector(0.,0.,-1.);
2297  }
2298  return snxt = 0.;
2299  }
2300  else if ( (eTheta > halfpi)
2301  && (t2 < 0.) && (p.z() <=0.) ) // leave
2302  {
2303  if( calcNorm ) { *validNorm = false; }
2304  return snxt = 0.;
2305  }
2306  }
2307  b = t2/t1;
2308  c = dist2ETheta/t1;
2309  d2 = b*b - c ;
2310  if ( (d2 <halfRmaxTolerance) && (d2 > -halfRmaxTolerance) )
2311  {
2312  d2 = 0.;
2313  }
2314  if ( d2 >= 0. )
2315  {
2316  d = std::sqrt(d2);
2317 
2318  if( eTheta < halfpi )
2319  {
2320  sd = -b - d; // First root
2321 
2322  if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.))
2323  || (sd < 0.) )
2324  {
2325  sd = -b + d ; // 2nd root
2326  }
2327  if( sd > halfRmaxTolerance )
2328  {
2329  if( sd < stheta )
2330  {
2331  stheta = sd;
2332  sidetheta = kETheta;
2333  }
2334  }
2335  }
2336  else // sTheta+fDTheta > pi/2, concave surface, no normal
2337  {
2338  sd = -b - d; // First root
2339 
2340  if ( ((std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.))
2341  || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > halfRmaxTolerance) ) )
2342  {
2343  sd = -b + d ; // 2nd root
2344  }
2345  if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= halfRmaxTolerance) )
2346  {
2347  if( sd < stheta )
2348  {
2349  stheta = sd;
2350  sidetheta = kETheta;
2351  }
2352  }
2353  }
2354  }
2355  }
2356  }
2357  }
2358 
2359  } // end theta intersections
2360 
2361  // Phi Intersection
2362 
2363  if ( !fFullPhiSphere )
2364  {
2365  if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
2366  {
2367  // pDist -ve when inside
2368 
2369  pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
2370  pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
2371 
2372  // Comp -ve when in direction of outwards normal
2373 
2374  compS = -sinSPhi*v.x()+cosSPhi*v.y() ;
2375  compE = sinEPhi*v.x()-cosEPhi*v.y() ;
2376  sidephi = kNull ;
2377 
2378  if ( (pDistS <= 0) && (pDistE <= 0) )
2379  {
2380  // Inside both phi *full* planes
2381 
2382  if ( compS < 0 )
2383  {
2384  sphi = pDistS/compS ;
2385  xi = p.x()+sphi*v.x() ;
2386  yi = p.y()+sphi*v.y() ;
2387 
2388  // Check intersection with correct half-plane (if not -> no intersect)
2389  //
2390  if( (std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance) )
2391  {
2392  vphi = std::atan2(v.y(),v.x());
2393  sidephi = kSPhi;
2394  if ( ( (fSPhi-halfAngTolerance) <= vphi)
2395  && ( (ePhi+halfAngTolerance) >= vphi) )
2396  {
2397  sphi = kInfinity;
2398  }
2399  }
2400  else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2401  {
2402  sphi=kInfinity;
2403  }
2404  else
2405  {
2406  sidephi = kSPhi ;
2407  if ( pDistS > -halfCarTolerance) { sphi = 0; } // Leave by sphi
2408  }
2409  }
2410  else { sphi = kInfinity; }
2411 
2412  if ( compE < 0 )
2413  {
2414  sphi2=pDistE/compE ;
2415  if (sphi2 < sphi) // Only check further if < starting phi intersection
2416  {
2417  xi = p.x()+sphi2*v.x() ;
2418  yi = p.y()+sphi2*v.y() ;
2419 
2420  // Check intersection with correct half-plane
2421  //
2422  if ((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
2423  {
2424  // Leaving via ending phi
2425  //
2426  vphi = std::atan2(v.y(),v.x()) ;
2427 
2428  if( !((fSPhi-halfAngTolerance <= vphi)
2429  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
2430  {
2431  sidephi = kEPhi;
2432  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
2433  else { sphi = 0.0; }
2434  }
2435  }
2436  else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
2437  {
2438  sidephi = kEPhi ;
2439  if ( pDistE <= -halfCarTolerance )
2440  {
2441  sphi=sphi2;
2442  }
2443  else
2444  {
2445  sphi = 0 ;
2446  }
2447  }
2448  }
2449  }
2450  }
2451  else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
2452  {
2453  if ( pDistS <= pDistE )
2454  {
2455  sidephi = kSPhi ;
2456  }
2457  else
2458  {
2459  sidephi = kEPhi ;
2460  }
2461  if ( fDPhi > pi )
2462  {
2463  if ( (compS < 0) && (compE < 0) ) { sphi = 0; }
2464  else { sphi = kInfinity; }
2465  }
2466  else
2467  {
2468  // if towards both >=0 then once inside (after error)
2469  // will remain inside
2470 
2471  if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
2472  else { sphi = 0; }
2473  }
2474  }
2475  else if ( (pDistS > 0) && (pDistE < 0) )
2476  {
2477  // Outside full starting plane, inside full ending plane
2478 
2479  if ( fDPhi > pi )
2480  {
2481  if ( compE < 0 )
2482  {
2483  sphi = pDistE/compE ;
2484  xi = p.x() + sphi*v.x() ;
2485  yi = p.y() + sphi*v.y() ;
2486 
2487  // Check intersection in correct half-plane
2488  // (if not -> not leaving phi extent)
2489  //
2490  if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2491  {
2492  vphi = std::atan2(v.y(),v.x());
2493  sidephi = kSPhi;
2494  if ( ( (fSPhi-halfAngTolerance) <= vphi)
2495  && ( (ePhi+halfAngTolerance) >= vphi) )
2496  {
2497  sphi = kInfinity;
2498  }
2499  }
2500  else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
2501  {
2502  sphi = kInfinity ;
2503  }
2504  else // Leaving via Ending phi
2505  {
2506  sidephi = kEPhi ;
2507  if ( pDistE > -halfCarTolerance ) { sphi = 0.; }
2508  }
2509  }
2510  else
2511  {
2512  sphi = kInfinity ;
2513  }
2514  }
2515  else
2516  {
2517  if ( compS >= 0 )
2518  {
2519  if ( compE < 0 )
2520  {
2521  sphi = pDistE/compE ;
2522  xi = p.x() + sphi*v.x() ;
2523  yi = p.y() + sphi*v.y() ;
2524 
2525  // Check intersection in correct half-plane
2526  // (if not -> remain in extent)
2527  //
2528  if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2529  {
2530  vphi = std::atan2(v.y(),v.x());
2531  sidephi = kSPhi;
2532  if ( ( (fSPhi-halfAngTolerance) <= vphi)
2533  && ( (ePhi+halfAngTolerance) >= vphi) )
2534  {
2535  sphi = kInfinity;
2536  }
2537  }
2538  else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
2539  {
2540  sphi=kInfinity;
2541  }
2542  else // otherwise leaving via Ending phi
2543  {
2544  sidephi = kEPhi ;
2545  }
2546  }
2547  else sphi=kInfinity;
2548  }
2549  else // leaving immediately by starting phi
2550  {
2551  sidephi = kSPhi ;
2552  sphi = 0 ;
2553  }
2554  }
2555  }
2556  else
2557  {
2558  // Must be pDistS < 0 && pDistE > 0
2559  // Inside full starting plane, outside full ending plane
2560 
2561  if ( fDPhi > pi )
2562  {
2563  if ( compS < 0 )
2564  {
2565  sphi=pDistS/compS;
2566  xi=p.x()+sphi*v.x();
2567  yi=p.y()+sphi*v.y();
2568 
2569  // Check intersection in correct half-plane
2570  // (if not -> not leaving phi extent)
2571  //
2572  if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
2573  {
2574  vphi = std::atan2(v.y(),v.x()) ;
2575  sidephi = kSPhi;
2576  if ( ( (fSPhi-halfAngTolerance) <= vphi)
2577  && ( (ePhi+halfAngTolerance) >= vphi) )
2578  {
2579  sphi = kInfinity;
2580  }
2581  }
2582  else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2583  {
2584  sphi = kInfinity ;
2585  }
2586  else // Leaving via Starting phi
2587  {
2588  sidephi = kSPhi ;
2589  if ( pDistS > -halfCarTolerance ) { sphi = 0; }
2590  }
2591  }
2592  else
2593  {
2594  sphi = kInfinity ;
2595  }
2596  }
2597  else
2598  {
2599  if ( compE >= 0 )
2600  {
2601  if ( compS < 0 )
2602  {
2603  sphi = pDistS/compS ;
2604  xi = p.x()+sphi*v.x() ;
2605  yi = p.y()+sphi*v.y() ;
2606 
2607  // Check intersection in correct half-plane
2608  // (if not -> remain in extent)
2609  //
2610  if((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
2611  {
2612  vphi = std::atan2(v.y(),v.x()) ;
2613  sidephi = kSPhi;
2614  if ( ( (fSPhi-halfAngTolerance) <= vphi)
2615  && ( (ePhi+halfAngTolerance) >= vphi) )
2616  {
2617  sphi = kInfinity;
2618  }
2619  }
2620  else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
2621  {
2622  sphi = kInfinity ;
2623  }
2624  else // otherwise leaving via Starting phi
2625  {
2626  sidephi = kSPhi ;
2627  }
2628  }
2629  else
2630  {
2631  sphi = kInfinity ;
2632  }
2633  }
2634  else // leaving immediately by ending
2635  {
2636  sidephi = kEPhi ;
2637  sphi = 0 ;
2638  }
2639  }
2640  }
2641  }
2642  else
2643  {
2644  // On z axis + travel not || to z axis -> if phi of vector direction
2645  // within phi of shape, Step limited by rmax, else Step =0
2646 
2647  if ( v.x() || v.y() )
2648  {
2649  vphi = std::atan2(v.y(),v.x()) ;
2650  if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
2651  {
2652  sphi = kInfinity;
2653  }
2654  else
2655  {
2656  sidephi = kSPhi ; // arbitrary
2657  sphi = 0 ;
2658  }
2659  }
2660  else // travel along z - no phi intersection
2661  {
2662  sphi = kInfinity ;
2663  }
2664  }
2665  if ( sphi < snxt ) // Order intersecttions
2666  {
2667  snxt = sphi ;
2668  side = sidephi ;
2669  }
2670  }
2671  if (stheta < snxt ) // Order intersections
2672  {
2673  snxt = stheta ;
2674  side = sidetheta ;
2675  }
2676 
2677  if (calcNorm) // Output switch operator
2678  {
2679  switch( side )
2680  {
2681  case kRMax:
2682  xi=p.x()+snxt*v.x();
2683  yi=p.y()+snxt*v.y();
2684  zi=p.z()+snxt*v.z();
2685  *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
2686  *validNorm=true;
2687  break;
2688 
2689  case kRMin:
2690  *validNorm=false; // Rmin is concave
2691  break;
2692 
2693  case kSPhi:
2694  if ( fDPhi <= pi ) // Normal to Phi-
2695  {
2696  *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
2697  *validNorm=true;
2698  }
2699  else { *validNorm=false; }
2700  break ;
2701 
2702  case kEPhi:
2703  if ( fDPhi <= pi ) // Normal to Phi+
2704  {
2705  *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
2706  *validNorm=true;
2707  }
2708  else { *validNorm=false; }
2709  break;
2710 
2711  case kSTheta:
2712  if( fSTheta == halfpi )
2713  {
2714  *n=G4ThreeVector(0.,0.,1.);
2715  *validNorm=true;
2716  }
2717  else if ( fSTheta > halfpi )
2718  {
2719  xi = p.x() + snxt*v.x();
2720  yi = p.y() + snxt*v.y();
2721  rho2=xi*xi+yi*yi;
2722  if (rho2)
2723  {
2724  rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
2725  *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
2726  -tanSTheta/std::sqrt(1+tanSTheta2));
2727  }
2728  else
2729  {
2730  *n = G4ThreeVector(0.,0.,1.);
2731  }
2732  *validNorm=true;
2733  }
2734  else { *validNorm=false; } // Concave STheta cone
2735  break;
2736 
2737  case kETheta:
2738  if( eTheta == halfpi )
2739  {
2740  *n = G4ThreeVector(0.,0.,-1.);
2741  *validNorm = true;
2742  }
2743  else if ( eTheta < halfpi )
2744  {
2745  xi=p.x()+snxt*v.x();
2746  yi=p.y()+snxt*v.y();
2747  rho2=xi*xi+yi*yi;
2748  if (rho2)
2749  {
2750  rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
2751  *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
2752  -tanETheta/std::sqrt(1+tanETheta2) );
2753  }
2754  else
2755  {
2756  *n = G4ThreeVector(0.,0.,-1.);
2757  }
2758  *validNorm=true;
2759  }
2760  else { *validNorm=false; } // Concave ETheta cone
2761  break;
2762 
2763  default:
2764  G4cout << G4endl;
2765  DumpInfo();
2766  std::ostringstream message;
2767  G4int oldprc = message.precision(16);
2768  message << "Undefined side for valid surface normal to solid."
2769  << G4endl
2770  << "Position:" << G4endl << G4endl
2771  << "p.x() = " << p.x()/mm << " mm" << G4endl
2772  << "p.y() = " << p.y()/mm << " mm" << G4endl
2773  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2774  << "Direction:" << G4endl << G4endl
2775  << "v.x() = " << v.x() << G4endl
2776  << "v.y() = " << v.y() << G4endl
2777  << "v.z() = " << v.z() << G4endl << G4endl
2778  << "Proposed distance :" << G4endl << G4endl
2779  << "snxt = " << snxt/mm << " mm" << G4endl;
2780  message.precision(oldprc);
2781  G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2782  "GeomSolids1002", JustWarning, message);
2783  break;
2784  }
2785  }
2786  if (snxt == kInfinity)
2787  {
2788  G4cout << G4endl;
2789  DumpInfo();
2790  std::ostringstream message;
2791  G4int oldprc = message.precision(16);
2792  message << "Logic error: snxt = kInfinity ???" << G4endl
2793  << "Position:" << G4endl << G4endl
2794  << "p.x() = " << p.x()/mm << " mm" << G4endl
2795  << "p.y() = " << p.y()/mm << " mm" << G4endl
2796  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2797  << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm
2798  << " mm" << G4endl << G4endl
2799  << "Direction:" << G4endl << G4endl
2800  << "v.x() = " << v.x() << G4endl
2801  << "v.y() = " << v.y() << G4endl
2802  << "v.z() = " << v.z() << G4endl << G4endl
2803  << "Proposed distance :" << G4endl << G4endl
2804  << "snxt = " << snxt/mm << " mm" << G4endl;
2805  message.precision(oldprc);
2806  G4Exception("G4Sphere::DistanceToOut(p,v,..)",
2807  "GeomSolids1002", JustWarning, message);
2808  }
2809 
2810  return snxt;
2811 }
2812 
2814 //
2815 // Calculate distance (<=actual) to closest surface of shape from inside
2816 
2817 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const
2818 {
2819  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
2820  G4double rho2,rds,rho;
2821  G4double pTheta,dTheta1 = kInfinity,dTheta2 = kInfinity;
2822  rho2=p.x()*p.x()+p.y()*p.y();
2823  rds=std::sqrt(rho2+p.z()*p.z());
2824  rho=std::sqrt(rho2);
2825 
2826 #ifdef G4CSGDEBUG
2827  if( Inside(p) == kOutside )
2828  {
2829  G4int old_prc = G4cout.precision(16);
2830  G4cout << G4endl;
2831  DumpInfo();
2832  G4cout << "Position:" << G4endl << G4endl ;
2833  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2834  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2835  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2836  G4cout.precision(old_prc) ;
2837  G4Exception("G4Sphere::DistanceToOut(p)",
2838  "GeomSolids1002", JustWarning, "Point p is outside !?" );
2839  }
2840 #endif
2841 
2842  // Distance to r shells
2843  //
2844  safeRMax = fRmax-rds;
2845  safe = safeRMax;
2846  if (fRmin)
2847  {
2848  safeRMin = rds-fRmin;
2849  safe = std::min( safeRMin, safeRMax );
2850  }
2851 
2852  // Distance to phi extent
2853  //
2854  if ( !fFullPhiSphere )
2855  {
2856  if (rho>0.0)
2857  {
2858  if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
2859  {
2860  safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
2861  }
2862  else
2863  {
2864  safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
2865  }
2866  }
2867  else
2868  {
2869  safePhi = 0.0; // Distance to both Phi surfaces (extended)
2870  }
2871  // Both cases above can be improved - in case fRMin > 0.0
2872  // although it may be costlier (good for precise, not fast version)
2873 
2874  safe= std::min(safe, safePhi);
2875  }
2876 
2877  // Distance to Theta extent
2878  //
2879  if ( !fFullThetaSphere )
2880  {
2881  if( rds > 0.0 )
2882  {
2883  pTheta=std::acos(p.z()/rds);
2884  if (pTheta<0) { pTheta+=pi; }
2885  if(fSTheta>0.)
2886  { dTheta1=pTheta-fSTheta;}
2887  if(eTheta<pi)
2888  { dTheta2=eTheta-pTheta;}
2889 
2890  safeTheta=rds*std::sin(std::min(dTheta1, dTheta2) );
2891  }
2892  else
2893  {
2894  safeTheta= 0.0;
2895  // An improvement will be to return negative answer if outside (TODO)
2896  }
2897  safe = std::min( safe, safeTheta );
2898  }
2899 
2900  if (safe<0.0) { safe=0; }
2901  // An improvement to return negative answer if outside (TODO)
2902 
2903  return safe;
2904 }
2905 
2907 //
2908 // Create a List containing the transformed vertices
2909 // Ordering [0-3] -fDz cross section
2910 // [4-7] +fDz cross section such that [0] is below [4],
2911 // [1] below [5] etc.
2912 // Note:
2913 // Caller has deletion resposibility
2914 // Potential improvement: For last slice, use actual ending angle
2915 // to avoid rounding error problems.
2916 
2917 G4ThreeVectorList*
2918 G4Sphere::CreateRotatedVertices( const G4AffineTransform& pTransform,
2919  G4int& noPolygonVertices ) const
2920 {
2921  G4ThreeVectorList *vertices;
2922  G4ThreeVector vertex;
2923  G4double meshAnglePhi,meshRMax,crossAnglePhi,
2924  coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
2925  G4double meshTheta,crossTheta,startTheta;
2926  G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
2927  G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;
2928 
2929  // Phi cross sections
2930 
2931  noPhiCrossSections = G4int(fDPhi/kMeshAngleDefault)+1;
2932 
2933  if (noPhiCrossSections<kMinMeshSections)
2934  {
2935  noPhiCrossSections=kMinMeshSections;
2936  }
2937  else if (noPhiCrossSections>kMaxMeshSections)
2938  {
2939  noPhiCrossSections=kMaxMeshSections;
2940  }
2941  meshAnglePhi=fDPhi/(noPhiCrossSections-1);
2942 
2943  // If complete in phi, set start angle such that mesh will be at fRMax
2944  // on the x axis. Will give better extent calculations when not rotated.
2945 
2946  if (fFullPhiSphere)
2947  {
2948  sAnglePhi = -meshAnglePhi*0.5;
2949  }
2950  else
2951  {
2952  sAnglePhi=fSPhi;
2953  }
2954 
2955  // Theta cross sections
2956 
2957  noThetaSections = G4int(fDTheta/kMeshAngleDefault)+1;
2958 
2959  if (noThetaSections<kMinMeshSections)
2960  {
2961  noThetaSections=kMinMeshSections;
2962  }
2963  else if (noThetaSections>kMaxMeshSections)
2964  {
2965  noThetaSections=kMaxMeshSections;
2966  }
2967  meshTheta=fDTheta/(noThetaSections-1);
2968 
2969  // If complete in Theta, set start angle such that mesh will be at fRMax
2970  // on the z axis. Will give better extent calculations when not rotated.
2971 
2972  if (fFullThetaSphere)
2973  {
2974  startTheta = -meshTheta*0.5;
2975  }
2976  else
2977  {
2978  startTheta=fSTheta;
2979  }
2980 
2981  meshRMax = (meshAnglePhi >= meshTheta) ?
2982  fRmax/std::cos(meshAnglePhi*0.5) : fRmax/std::cos(meshTheta*0.5);
2983  G4double* cosCrossTheta = new G4double[noThetaSections];
2984  G4double* sinCrossTheta = new G4double[noThetaSections];
2985  vertices=new G4ThreeVectorList();
2986  if (vertices && cosCrossTheta && sinCrossTheta)
2987  {
2988  vertices->reserve(noPhiCrossSections*(noThetaSections*2));
2989  for (crossSectionPhi=0;
2990  crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
2991  {
2992  crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
2993  coscrossAnglePhi=std::cos(crossAnglePhi);
2994  sincrossAnglePhi=std::sin(crossAnglePhi);
2995  for (crossSectionTheta=0;
2996  crossSectionTheta<noThetaSections;crossSectionTheta++)
2997  {
2998  // Compute coordinates of cross section at section crossSectionPhi
2999  //
3000  crossTheta=startTheta+crossSectionTheta*meshTheta;
3001  cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
3002  sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
3003 
3004  rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
3005  rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
3006  rMinZ=fRmin*cosCrossTheta[crossSectionTheta];
3007 
3008  vertex=G4ThreeVector(rMinX,rMinY,rMinZ);
3009  vertices->push_back(pTransform.TransformPoint(vertex));
3010 
3011  } // Theta forward
3012 
3013  for (crossSectionTheta=noThetaSections-1;
3014  crossSectionTheta>=0; crossSectionTheta--)
3015  {
3016  rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
3017  rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
3018  rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
3019 
3020  vertex=G4ThreeVector(rMaxX,rMaxY,rMaxZ);
3021  vertices->push_back(pTransform.TransformPoint(vertex));
3022 
3023  } // Theta back
3024  } // Phi
3025  noPolygonVertices = noThetaSections*2 ;
3026  }
3027  else
3028  {
3029  DumpInfo();
3030  G4Exception("G4Sphere::CreateRotatedVertices()",
3031  "GeomSolids0003", FatalException,
3032  "Error in allocation of vertices. Out of memory !");
3033  }
3034 
3035  delete [] cosCrossTheta;
3036  delete [] sinCrossTheta;
3037 
3038  return vertices;
3039 }
3040 
3042 //
3043 // G4EntityType
3044 
3045 G4GeometryType G4Sphere::GetEntityType() const
3046 {
3047  return G4String("G4Sphere");
3048 }
3049 
3051 //
3052 // Make a clone of the object
3053 //
3054 G4VSolid* G4Sphere::Clone() const
3055 {
3056  return new G4Sphere(*this);
3057 }
3058 
3060 //
3061 // Stream object contents to an output stream
3062 
3063 std::ostream& G4Sphere::StreamInfo( std::ostream& os ) const
3064 {
3065  G4int oldprc = os.precision(16);
3066  os << "-----------------------------------------------------------\n"
3067  << " *** Dump for solid - " << GetName() << " ***\n"
3068  << " ===================================================\n"
3069  << " Solid type: G4Sphere\n"
3070  << " Parameters: \n"
3071  << " inner radius: " << fRmin/mm << " mm \n"
3072  << " outer radius: " << fRmax/mm << " mm \n"
3073  << " starting phi of segment : " << fSPhi/degree << " degrees \n"
3074  << " delta phi of segment : " << fDPhi/degree << " degrees \n"
3075  << " starting theta of segment: " << fSTheta/degree << " degrees \n"
3076  << " delta theta of segment : " << fDTheta/degree << " degrees \n"
3077  << "-----------------------------------------------------------\n";
3078  os.precision(oldprc);
3079 
3080  return os;
3081 }
3082 
3084 //
3085 // GetPointOnSurface
3086 
3087 G4ThreeVector G4Sphere::GetPointOnSurface() const
3088 {
3089  G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
3090  G4double height1, height2, slant1, slant2, costheta, sintheta, rRand;
3091 
3092  height1 = (fRmax-fRmin)*cosSTheta;
3093  height2 = (fRmax-fRmin)*cosETheta;
3094  slant1 = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
3095  slant2 = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
3096  rRand = GetRadiusInRing(fRmin,fRmax);
3097 
3098  aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
3099  aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
3100  aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
3101  aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
3102  aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
3103 
3104  phi = G4RandFlat::shoot(fSPhi, ePhi);
3105  cosphi = std::cos(phi);
3106  sinphi = std::sin(phi);
3107  costheta = G4RandFlat::shoot(cosETheta,cosSTheta);
3108  sintheta = std::sqrt(1.-sqr(costheta));
3109 
3110  if(fFullPhiSphere) { aFiv = 0; }
3111  if(fSTheta == 0) { aThr=0; }
3112  if(eTheta == pi) { aFou = 0; }
3113  if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
3114  if(eTheta == halfpi) { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }
3115 
3116  chose = G4RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
3117  if( (chose>=0.) && (chose<aOne) )
3118  {
3119  return G4ThreeVector(fRmax*sintheta*cosphi,
3120  fRmax*sintheta*sinphi, fRmax*costheta);
3121  }
3122  else if( (chose>=aOne) && (chose<aOne+aTwo) )
3123  {
3124  return G4ThreeVector(fRmin*sintheta*cosphi,
3125  fRmin*sintheta*sinphi, fRmin*costheta);
3126  }
3127  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
3128  {
3129  if (fSTheta != halfpi)
3130  {
3131  zRand = G4RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
3132  return G4ThreeVector(tanSTheta*zRand*cosphi,
3133  tanSTheta*zRand*sinphi,zRand);
3134  }
3135  else
3136  {
3137  return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
3138  }
3139  }
3140  else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
3141  {
3142  if(eTheta != halfpi)
3143  {
3144  zRand = G4RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
3145  return G4ThreeVector (tanETheta*zRand*cosphi,
3146  tanETheta*zRand*sinphi,zRand);
3147  }
3148  else
3149  {
3150  return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
3151  }
3152  }
3153  else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
3154  {
3155  return G4ThreeVector(rRand*sintheta*cosSPhi,
3156  rRand*sintheta*sinSPhi,rRand*costheta);
3157  }
3158  else
3159  {
3160  return G4ThreeVector(rRand*sintheta*cosEPhi,
3161  rRand*sintheta*sinEPhi,rRand*costheta);
3162  }
3163 }
3164 
3166 //
3167 // GetSurfaceArea
3168 
3169 G4double G4Sphere::GetSurfaceArea()
3170 {
3171  if(fSurfaceArea != 0.) {;}
3172  else
3173  {
3174  G4double Rsq=fRmax*fRmax;
3175  G4double rsq=fRmin*fRmin;
3176 
3177  fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
3178  if(!fFullPhiSphere)
3179  {
3180  fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
3181  }
3182  if(fSTheta >0)
3183  {
3184  G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
3185  + std::pow(cosSTheta,2));
3186  if(fDPhi>pi)
3187  {
3188  fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
3189  }
3190  else
3191  {
3192  fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
3193  }
3194  }
3195  if(eTheta < pi)
3196  {
3197  G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
3198  + std::pow(cosETheta,2));
3199  if(fDPhi>pi)
3200  {
3201  fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
3202  }
3203  else
3204  {
3205  fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
3206  }
3207  }
3208  }
3209  return fSurfaceArea;
3210 }
3211 
3213 //
3214 // Methods for visualisation
3215 
3216 G4VisExtent G4Sphere::GetExtent() const
3217 {
3218  return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
3219 }
3220 
3221 
3222 void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
3223 {
3224  scene.AddSolid (*this);
3225 }
3226 
3227 G4Polyhedron* G4Sphere::CreatePolyhedron () const
3228 {
3229  return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
3230 }
3231 
3232 #endif
kNRMin
Definition: G4Sphere.cc:82
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
t1
TTree * t1
Definition: plottest35.C:26
G4Sphere::fDTheta
G4double fDTheta
Definition: G4Sphere.hh:234
G4Sphere::fDPhi
G4double fDPhi
Definition: G4Sphere.hh:234
G4Sphere::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Sphere.cc:888
G4Sphere::halfCarTolerance
G4double halfCarTolerance
Definition: G4Sphere.hh:252
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4Sphere::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Sphere.cc:472
G4Sphere::fEpsilon
G4double fEpsilon
Definition: G4Sphere.hh:229
halfpi
static const double halfpi
Definition: G4SIunits.hh:76
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VGraphicsScene.hh
d
Float_t d
Definition: advanced/microbeam/plot.C:237
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Sphere::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Sphere.cc:3087
G4Sphere::tanETheta2
G4double tanETheta2
Definition: G4Sphere.hh:243
G4VoxelLimits.hh
G4Sphere::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4Sphere.cc:3045
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:111
in
ifstream in
Definition: comparison.C:7
kNSTheta
Definition: G4Sphere.cc:82
G4Sphere::fSTheta
G4double fSTheta
Definition: G4Sphere.hh:234
kRMax
Definition: G4Sphere.cc:78
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
norm
Float_t norm
Definition: comparison_ascii.C:71
G4Sphere::fRmin
G4double fRmin
Definition: G4Sphere.hh:234
G4Sphere::fRminTolerance
G4double fRminTolerance
Definition: G4Sphere.hh:229
G4Sphere::hDPhi
G4double hDPhi
Definition: G4Sphere.hh:238
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
kEPhi
Definition: G4Sphere.cc:78
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
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4Sphere::G4Sphere
G4Sphere(const G4String &pName, G4double pRmin, G4double pRmax, G4double pSPhi, G4double pDPhi, G4double pSTheta, G4double pDTheta)
Definition: G4Sphere.cc:89
G4VSolid::GetName
G4String GetName() const
G4VSolid::DumpInfo
void DumpInfo() const
G4Sphere::cosSTheta
G4double cosSTheta
Definition: G4Sphere.hh:243
ENorm
ENorm
Definition: G4Cons.cc:75
G4Sphere::sinSPhi
G4double sinSPhi
Definition: G4Sphere.hh:238
s
static const double s
Definition: G4SIunits.hh:168
G4Sphere::~G4Sphere
~G4Sphere()
Definition: G4Sphere.cc:144
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
G4Sphere::cosHDPhiOT
G4double cosHDPhiOT
Definition: G4Sphere.hh:238
G4VisExtent.hh
n
Char_t n[5]
Definition: comparison_ascii.C:55
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4CSGSolid::fSurfaceArea
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4cout
G4GLOB_DLL std::ostream G4cout
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
G4Sphere::sinEPhi
G4double sinEPhi
Definition: G4Sphere.hh:238
kRMin
Definition: G4Sphere.cc:78
G4Sphere::CheckThetaAngles
void CheckThetaAngles(G4double sTheta, G4double dTheta)
G4Sphere::fRmax
G4double fRmax
Definition: G4Sphere.hh:234
G4bool
bool G4bool
Definition: G4Types.hh:79
G4Sphere::tanSTheta
G4double tanSTheta
Definition: G4Sphere.hh:243
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
twopi
static const double twopi
Definition: G4SIunits.hh:75
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
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition: G4GeometryTolerance.cc:92
Randomize.hh
G4VoxelLimits::GetMinExtent
G4double GetMinExtent(const EAxis pAxis) const
kSTheta
Definition: G4Sphere.cc:78
G4Sphere::sinCPhi
G4double sinCPhi
Definition: G4Sphere.hh:238
G4Sphere::CreateRotatedVertices
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform, G4int &noPolygonVertices) const
Definition: G4Sphere.cc:2918
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
CLHEP::Hep3Vector::x
double x() const
G4Sphere::operator=
G4Sphere & operator=(const G4Sphere &rhs)
Definition: G4Sphere.cc:178
G4Sphere::fRmaxTolerance
G4double fRmaxTolerance
Definition: G4Sphere.hh:229
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
G4Sphere::tanSTheta2
G4double tanSTheta2
Definition: G4Sphere.hh:243
G4Sphere::kRadTolerance
G4double kRadTolerance
Definition: G4Sphere.hh:229
G4AffineTransform.hh
pi
static const double pi
Definition: G4SIunits.hh:74
kNETheta
Definition: G4Sphere.cc:82
G4Sphere::tanETheta
G4double tanETheta
Definition: G4Sphere.hh:243
CLHEP::Hep3Vector::y
double y() const
G4Sphere::eTheta
G4double eTheta
Definition: G4Sphere.hh:243
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::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Sphere.cc:586
G4Sphere::sinETheta
G4double sinETheta
Definition: G4Sphere.hh:243
G4Sphere::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Sphere.cc:3063
G4Sphere::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Sphere.cc:216
G4Sphere::cPhi
G4double cPhi
Definition: G4Sphere.hh:238
CLHEP::Hep3Vector::z
double z() const
CLHEP
Definition: DoubConv.h:18
kNSPhi
Definition: G4Sphere.cc:82
FatalException
Definition: G4ExceptionSeverity.hh:60
degree
static const double degree
Definition: G4SIunits.hh:143
G4Sphere::GetSurfaceArea
G4double GetSurfaceArea()
Definition: G4Sphere.cc:3169
G4Sphere::kAngTolerance
G4double kAngTolerance
Definition: G4Sphere.hh:229
G4VoxelLimits::GetMaxExtent
G4double GetMaxExtent(const EAxis pAxis) const
t2
TTree * t2
Definition: plottest35.C:36
kNull
Definition: G4Sphere.cc:78
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4Sphere::cosEPhi
G4double cosEPhi
Definition: G4Sphere.hh:238
ESide
ESide
Definition: G4Cons.cc:71
sqr
T sqr(const T &x)
Definition: templates.hh:145
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:1932
G4Sphere::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Sphere.cc:3222
G4double
double G4double
Definition: G4Types.hh:76
G4Sphere::sinSTheta
G4double sinSTheta
Definition: G4Sphere.hh:243
G4VSolid::CalculateClippedPolygonExtent
void CalculateClippedPolygonExtent(G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:427
d2
static const G4double d2
Definition: G4KleinNishinaCompton.cc:66
kNEPhi
Definition: G4Sphere.cc:82
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4GeometryTolerance.hh
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
G4Sphere::GetExtent
G4VisExtent GetExtent() const
Definition: G4Sphere.cc:3216
G4Sphere.hh
mm
static const double mm
Definition: G4SIunits.hh:114
G4Sphere::halfAngTolerance
G4double halfAngTolerance
Definition: G4Sphere.hh:252
G4Sphere::cosHDPhiIT
G4double cosHDPhiIT
Definition: G4Sphere.hh:238
kMeshAngleDefault
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42
G4Sphere::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Sphere.cc:715
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
G4Sphere::fFullThetaSphere
G4bool fFullThetaSphere
Definition: G4Sphere.hh:248
kMaxMeshSections
const G4int kMaxMeshSections
Definition: meshdefs.hh:46
G4Sphere::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Sphere.cc:3227
G4Sphere::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Sphere.cc:227
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:95
G4Sphere::cosCPhi
G4double cosCPhi
Definition: G4Sphere.hh:238
G4Sphere::Clone
G4VSolid * Clone() const
Definition: G4Sphere.cc:3054