Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: