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

Last change on this file since 833 was 831, checked in by garnier, 17 years ago

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