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