source: trunk/source/geometry/solids/CSG/src/G4Sphere.cc@ 1316

Last change on this file since 1316 was 1228, checked in by garnier, 16 years ago

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