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