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