source: trunk/source/geometry/solids/CSG/src/G4Torus.cc@ 1325

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

update geant4.9.3 tag
File size: 51.9 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: G4Torus.cc,v 1.65 2009/11/26 10:31:06 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30//
31// class G4Torus
32//
33// Implementation
34//
35// 02.10.07 T.Nikitina: Bug fixed in SolveNumericJT(), b.969:segmentation fault.
36// rootsrefined is used only if the number of refined roots
37// is the same as for primary roots.
38// 02.10.07 T.Nikitina: Bug fixed in CalculateExtent() for case of non-rotated
39// full-phi torus:protect against negative value for sqrt,
40// correct formula for delta.
41// 20.11.05 V.Grichine: Bug fixed in Inside(p) for phi sections, b.810
42// 25.08.05 O.Link: new methods for DistanceToIn/Out using JTPolynomialSolver
43// 07.06.05 V.Grichine: SurfaceNormal(p) for rho=0, Constructor as G4Cons
44// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
45// 18.03.04 V.Grichine: bug fixed in DistanceToIn(p)
46// 11.01.01 E.Medernach: Use G4PolynomialSolver to find roots
47// 03.10.00 E.Medernach: SafeNewton added
48// 31.08.00 E.Medernach: numerical computation of roots wuth bounding
49// volume technique
50// 26.05.00 V.Grichine: new fuctions developed by O.Cremonesi were added
51// 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
52// 19.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
53// 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...)
54// 30.10.96 V.Grichine: first implementation with G4Tubs elements in Fs
55//
56
57#include "G4Torus.hh"
58
59#include "G4VoxelLimits.hh"
60#include "G4AffineTransform.hh"
61#include "G4GeometryTolerance.hh"
62#include "G4JTPolynomialSolver.hh"
63
64#include "G4VPVParameterisation.hh"
65
66#include "meshdefs.hh"
67
68#include "Randomize.hh"
69
70#include "G4VGraphicsScene.hh"
71#include "G4Polyhedron.hh"
72#include "G4NURBS.hh"
73#include "G4NURBStube.hh"
74#include "G4NURBScylinder.hh"
75#include "G4NURBStubesector.hh"
76
77using namespace CLHEP;
78
79///////////////////////////////////////////////////////////////
80//
81// Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
82// - note if pdphi>2PI then reset to 2PI
83
84G4Torus::G4Torus( const G4String &pName,
85 G4double pRmin,
86 G4double pRmax,
87 G4double pRtor,
88 G4double pSPhi,
89 G4double pDPhi)
90 : G4CSGSolid(pName)
91{
92 SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi);
93}
94
95////////////////////////////////////////////////////////////////////////////
96//
97//
98
99void
100G4Torus::SetAllParameters( G4double pRmin,
101 G4double pRmax,
102 G4double pRtor,
103 G4double pSPhi,
104 G4double pDPhi )
105{
106 fCubicVolume = 0.;
107 fSurfaceArea = 0.;
108 fpPolyhedron = 0;
109
110 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
111 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
112
113 if ( pRtor >= pRmax+1.e3*kCarTolerance ) // Check swept radius, as in G4Cons
114 {
115 fRtor = pRtor ;
116 }
117 else
118 {
119 G4cerr << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl
120 << " Invalid swept radius !" << G4endl
121 << "pRtor = " << pRtor << ", pRmax = " << pRmax << G4endl;
122 G4Exception("G4Torus::SetAllParameters()",
123 "InvalidSetup", FatalException, "Invalid swept radius.");
124 }
125
126 // Check radii, as in G4Cons
127 //
128 if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
129 {
130 if (pRmin >= 1.e2*kCarTolerance) { fRmin = pRmin ; }
131 else { fRmin = 0.0 ; }
132 fRmax = pRmax ;
133 }
134 else
135 {
136 G4cerr << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl
137 << " Invalid values for radii !" << G4endl
138 << " pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl;
139 G4Exception("G4Torus::SetAllParameters()",
140 "InvalidSetup", FatalException, "Invalid radii.");
141 }
142
143 // Check angles
144 //
145 if ( pDPhi >= twopi ) { fDPhi = twopi ; }
146 else
147 {
148 if (pDPhi > 0) { fDPhi = pDPhi ; }
149 else
150 {
151 G4cerr << "ERROR - G4Torus::SetAllParameters(): " << GetName() << G4endl
152 << " Negative Z delta-Phi ! - "
153 << pDPhi << G4endl;
154 G4Exception("G4Torus::SetAllParameters()",
155 "InvalidSetup", FatalException, "Invalid dphi.");
156 }
157 }
158
159 // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0
160 //
161 fSPhi = pSPhi;
162
163 if (fSPhi < 0) { fSPhi = twopi-std::fmod(std::fabs(fSPhi),twopi) ; }
164 else { fSPhi = std::fmod(fSPhi,twopi) ; }
165
166 if (fSPhi+fDPhi > twopi) { fSPhi-=twopi ; }
167}
168
169///////////////////////////////////////////////////////////////////////
170//
171// Fake default constructor - sets only member data and allocates memory
172// for usage restricted to object persistency.
173//
174G4Torus::G4Torus( __void__& a )
175 : G4CSGSolid(a)
176{
177}
178
179//////////////////////////////////////////////////////////////////////
180//
181// Destructor
182
183G4Torus::~G4Torus()
184{}
185
186//////////////////////////////////////////////////////////////////////
187//
188// Dispatch to parameterisation for replication mechanism dimension
189// computation & modification.
190
191void G4Torus::ComputeDimensions( G4VPVParameterisation* p,
192 const G4int n,
193 const G4VPhysicalVolume* pRep )
194{
195 p->ComputeDimensions(*this,n,pRep);
196}
197
198
199
200////////////////////////////////////////////////////////////////////////////////
201//
202// Calculate the real roots to torus surface.
203// Returns negative solutions as well.
204
205std::vector<G4double> G4Torus::TorusRootsJT( const G4ThreeVector& p,
206 const G4ThreeVector& v,
207 G4double r ) const
208{
209
210 G4int i, num ;
211 G4double c[5], sr[4], si[4] ;
212 std::vector<G4double> roots ;
213
214 G4double Rtor2 = fRtor*fRtor, r2 = r*r ;
215
216 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
217 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
218
219 c[0] = 1.0 ;
220 c[1] = 4*pDotV ;
221 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - r2 + 2*Rtor2*v.z()*v.z()) ;
222 c[3] = 4*(pDotV*(pRad2 - Rtor2 - r2) + 2*Rtor2*p.z()*v.z()) ;
223 c[4] = pRad2*pRad2 - 2*pRad2*(Rtor2+r2)
224 + 4*Rtor2*p.z()*p.z() + (Rtor2-r2)*(Rtor2-r2) ;
225
226 G4JTPolynomialSolver torusEq;
227
228 num = torusEq.FindRoots( c, 4, sr, si );
229
230 for ( i = 0; i < num; i++ )
231 {
232 if( si[i] == 0. ) { roots.push_back(sr[i]) ; } // store real roots
233 }
234
235 std::sort(roots.begin() , roots.end() ) ; // sorting with <
236
237 return roots;
238}
239
240//////////////////////////////////////////////////////////////////////////////
241//
242// Interface for DistanceToIn and DistanceToOut.
243// Calls TorusRootsJT and returns the smalles possible distance to
244// the surface.
245// Attention: Difference in DistanceToIn/Out for points p on the surface.
246
247G4double G4Torus::SolveNumericJT( const G4ThreeVector& p,
248 const G4ThreeVector& v,
249 G4double r,
250 G4bool IsDistanceToIn ) const
251{
252 G4double bigdist = 10*mm ;
253 G4double tmin = kInfinity ;
254 G4double t, scal ;
255
256 // calculate the distances to the intersections with the Torus
257 // from a given point p and direction v.
258 //
259 std::vector<G4double> roots ;
260 std::vector<G4double> rootsrefined ;
261 roots = TorusRootsJT(p,v,r) ;
262
263 G4ThreeVector ptmp ;
264
265 // determine the smallest non-negative solution
266 //
267 for ( size_t k = 0 ; k<roots.size() ; k++ )
268 {
269 t = roots[k] ;
270
271 if ( t < -0.5*kCarTolerance ) { continue ; } // skip negative roots
272
273 if ( t > bigdist && t<kInfinity ) // problem with big distances
274 {
275 ptmp = p + t*v ;
276 rootsrefined = TorusRootsJT(ptmp,v,r) ;
277 if ( rootsrefined.size()==roots.size() )
278 {
279 t = t + rootsrefined[k] ;
280 }
281 }
282
283 ptmp = p + t*v ; // calculate the position of the proposed intersection
284
285 G4double theta = std::atan2(ptmp.y(),ptmp.x());
286
287 if ( fSPhi >= 0 )
288 {
289 if ( theta < - kAngTolerance*0.5 ) { theta += twopi; }
290 if ( (std::abs(theta) < kAngTolerance*0.5)
291 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
292 {
293 theta += twopi ; // 0 <= theta < 2pi
294 }
295 }
296 if ((fSPhi <= -pi )&&(theta>kAngTolerance*0.5)) { theta = theta-twopi; }
297
298 // We have to verify if this root is inside the region between
299 // fSPhi and fSPhi + fDPhi
300 //
301 if ( (theta - fSPhi >= - kAngTolerance*0.5)
302 && (theta - (fSPhi + fDPhi) <= kAngTolerance*0.5) )
303 {
304 // check if P is on the surface, and called from DistanceToIn
305 // DistanceToIn has to return 0.0 if particle is going inside the solid
306
307 if ( IsDistanceToIn == true )
308 {
309 if (std::fabs(t) < 0.5*kCarTolerance )
310 {
311 // compute scalar product at position p : v.n
312 // ( n taken from SurfaceNormal, not normalized )
313
314 scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
315 + p.y()*p.y())),
316 p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
317 + p.y()*p.y())),
318 p.z() );
319
320 // change sign in case of inner radius
321 //
322 if ( r == GetRmin() ) { scal = -scal ; }
323 if ( scal < 0 ) { return 0.0 ; }
324 }
325 }
326
327 // check if P is on the surface, and called from DistanceToOut
328 // DistanceToIn has to return 0.0 if particle is leaving the solid
329
330 if ( IsDistanceToIn == false )
331 {
332 if (std::fabs(t) < 0.5*kCarTolerance )
333 {
334 // compute scalar product at position p : v.n
335 //
336 scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
337 + p.y()*p.y())),
338 p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
339 + p.y()*p.y())),
340 p.z() );
341
342 // change sign in case of inner radius
343 //
344 if ( r == GetRmin() ) { scal = -scal ; }
345 if ( scal > 0 ) { return 0.0 ; }
346 }
347 }
348
349 // check if distance is larger than 1/2 kCarTolerance
350 //
351 if( t > 0.5*kCarTolerance )
352 {
353 tmin = t ;
354 return tmin ;
355 }
356 }
357 }
358
359 return tmin;
360}
361
362/////////////////////////////////////////////////////////////////////////////
363//
364// Calculate extent under transform and specified limit
365
366G4bool G4Torus::CalculateExtent( const EAxis pAxis,
367 const G4VoxelLimits& pVoxelLimit,
368 const G4AffineTransform& pTransform,
369 G4double& pMin, G4double& pMax) const
370{
371 if ((!pTransform.IsRotated()) && (fDPhi==twopi) && (fRmin==0))
372 {
373 // Special case handling for unrotated solid torus
374 // Compute x/y/z mins and maxs for bounding box respecting limits,
375 // with early returns if outside limits. Then switch() on pAxis,
376 // and compute exact x and y limit for x/y case
377
378 G4double xoffset,xMin,xMax;
379 G4double yoffset,yMin,yMax;
380 G4double zoffset,zMin,zMax;
381
382 G4double RTorus,delta,diff1,diff2,maxDiff,newMin,newMax;
383 G4double xoff1,xoff2,yoff1,yoff2;
384
385 xoffset = pTransform.NetTranslation().x();
386 xMin = xoffset - fRmax - fRtor ;
387 xMax = xoffset + fRmax + fRtor ;
388
389 if (pVoxelLimit.IsXLimited())
390 {
391 if ( (xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance)
392 || (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance) )
393 return false ;
394 else
395 {
396 if (xMin < pVoxelLimit.GetMinXExtent())
397 {
398 xMin = pVoxelLimit.GetMinXExtent() ;
399 }
400 if (xMax > pVoxelLimit.GetMaxXExtent())
401 {
402 xMax = pVoxelLimit.GetMaxXExtent() ;
403 }
404 }
405 }
406 yoffset = pTransform.NetTranslation().y();
407 yMin = yoffset - fRmax - fRtor ;
408 yMax = yoffset + fRmax + fRtor ;
409
410 if (pVoxelLimit.IsYLimited())
411 {
412 if ( (yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance)
413 || (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance) )
414 {
415 return false ;
416 }
417 else
418 {
419 if (yMin < pVoxelLimit.GetMinYExtent() )
420 {
421 yMin = pVoxelLimit.GetMinYExtent() ;
422 }
423 if (yMax > pVoxelLimit.GetMaxYExtent() )
424 {
425 yMax = pVoxelLimit.GetMaxYExtent() ;
426 }
427 }
428 }
429 zoffset = pTransform.NetTranslation().z() ;
430 zMin = zoffset - fRmax ;
431 zMax = zoffset + fRmax ;
432
433 if (pVoxelLimit.IsZLimited())
434 {
435 if ( (zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance)
436 || (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance) )
437 {
438 return false ;
439 }
440 else
441 {
442 if (zMin < pVoxelLimit.GetMinZExtent() )
443 {
444 zMin = pVoxelLimit.GetMinZExtent() ;
445 }
446 if (zMax > pVoxelLimit.GetMaxZExtent() )
447 {
448 zMax = pVoxelLimit.GetMaxZExtent() ;
449 }
450 }
451 }
452
453 // Known to cut cylinder
454
455 switch (pAxis)
456 {
457 case kXAxis:
458 yoff1=yoffset-yMin;
459 yoff2=yMax-yoffset;
460 if ( yoff1 >= 0 && yoff2 >= 0 )
461 {
462 // Y limits cross max/min x => no change
463 //
464 pMin = xMin ;
465 pMax = xMax ;
466 }
467 else
468 {
469 // Y limits don't cross max/min x => compute max delta x,
470 // hence new mins/maxs
471 //
472 RTorus=fRmax+fRtor;
473 delta = RTorus*RTorus - yoff1*yoff1;
474 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
475 delta = RTorus*RTorus - yoff2*yoff2;
476 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
477 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
478 newMin = xoffset - maxDiff ;
479 newMax = xoffset + maxDiff ;
480 pMin = (newMin < xMin) ? xMin : newMin ;
481 pMax = (newMax > xMax) ? xMax : newMax ;
482 }
483 break;
484
485 case kYAxis:
486 xoff1 = xoffset - xMin ;
487 xoff2 = xMax - xoffset ;
488 if (xoff1 >= 0 && xoff2 >= 0 )
489 {
490 // X limits cross max/min y => no change
491 //
492 pMin = yMin ;
493 pMax = yMax ;
494 }
495 else
496 {
497 // X limits don't cross max/min y => compute max delta y,
498 // hence new mins/maxs
499 //
500 RTorus=fRmax+fRtor;
501 delta = RTorus*RTorus - xoff1*xoff1;
502 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
503 delta = RTorus*RTorus - xoff2*xoff2;
504 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
505 maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
506 newMin = yoffset - maxDiff ;
507 newMax = yoffset + maxDiff ;
508 pMin = (newMin < yMin) ? yMin : newMin ;
509 pMax = (newMax > yMax) ? yMax : newMax ;
510 }
511 break;
512
513 case kZAxis:
514 pMin=zMin;
515 pMax=zMax;
516 break;
517
518 default:
519 break;
520 }
521 pMin -= kCarTolerance ;
522 pMax += kCarTolerance ;
523
524 return true;
525 }
526 else
527 {
528 G4int i, noEntries, noBetweenSections4 ;
529 G4bool existsAfterClip = false ;
530
531 // Calculate rotated vertex coordinates
532
533 G4ThreeVectorList *vertices ;
534 G4int noPolygonVertices ; // will be 4
535 vertices = CreateRotatedVertices(pTransform,noPolygonVertices) ;
536
537 pMin = +kInfinity ;
538 pMax = -kInfinity ;
539
540 noEntries = vertices->size() ;
541 noBetweenSections4 = noEntries - noPolygonVertices ;
542
543 for (i=0;i<noEntries;i+=noPolygonVertices)
544 {
545 ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
546 }
547 for (i=0;i<noBetweenSections4;i+=noPolygonVertices)
548 {
549 ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
550 }
551 if (pMin!=kInfinity||pMax!=-kInfinity)
552 {
553 existsAfterClip = true ; // Add 2*tolerance to avoid precision troubles
554 pMin -= kCarTolerance ;
555 pMax += kCarTolerance ;
556 }
557 else
558 {
559 // Check for case where completely enveloping clipping volume
560 // If point inside then we are confident that the solid completely
561 // envelopes the clipping volume. Hence set min/max extents according
562 // to clipping volume extents along the specified axis.
563
564 G4ThreeVector clipCentre(
565 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
566 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
567 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ;
568
569 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
570 {
571 existsAfterClip = true ;
572 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
573 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
574 }
575 }
576 delete vertices;
577 return existsAfterClip;
578 }
579}
580
581//////////////////////////////////////////////////////////////////////////////
582//
583// Return whether point inside/outside/on surface
584
585EInside G4Torus::Inside( const G4ThreeVector& p ) const
586{
587 G4double r2, pt2, pPhi, tolRMin, tolRMax ;
588
589 EInside in = kOutside ;
590 // General precals
591 r2 = p.x()*p.x() + p.y()*p.y() ;
592 pt2 = r2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*std::sqrt(r2) ;
593
594 if (fRmin) tolRMin = fRmin + kRadTolerance*0.5 ;
595 else tolRMin = 0 ;
596
597 tolRMax = fRmax - kRadTolerance*0.5;
598
599 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
600 {
601 if ( fDPhi == twopi || pt2 == 0 ) // on torus swept axis
602 {
603 in = kInside ;
604 }
605 else
606 {
607 // Try inner tolerant phi boundaries (=>inside)
608 // if not inside, try outer tolerant phi boundaries
609
610 pPhi = std::atan2(p.y(),p.x()) ;
611
612 if ( pPhi < -kAngTolerance*0.5 ) { pPhi += twopi ; } // 0<=pPhi<2pi
613 if ( fSPhi >= 0 )
614 {
615 if ( (std::abs(pPhi) < kAngTolerance*0.5)
616 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
617 {
618 pPhi += twopi ; // 0 <= pPhi < 2pi
619 }
620 if ( (pPhi >= fSPhi + kAngTolerance*0.5)
621 && (pPhi <= fSPhi + fDPhi - kAngTolerance*0.5) )
622 {
623 in = kInside ;
624 }
625 else if ( (pPhi >= fSPhi - kAngTolerance*0.5)
626 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
627 {
628 in = kSurface ;
629 }
630 }
631 else // fSPhi < 0
632 {
633 if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
634 && (pPhi >= fSPhi + fDPhi + kAngTolerance*0.5) ) {;}
635 else
636 {
637 in = kSurface ;
638 }
639 }
640 }
641 }
642 else // Try generous boundaries
643 {
644 tolRMin = fRmin - kRadTolerance*0.5 ;
645 tolRMax = fRmax + kRadTolerance*0.5 ;
646
647 if (tolRMin < 0 ) { tolRMin = 0 ; }
648
649 if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
650 {
651 if ( (fDPhi == twopi) || (pt2 == 0) ) // Continuous in phi or on z-axis
652 {
653 in = kSurface ;
654 }
655 else // Try outer tolerant phi boundaries only
656 {
657 pPhi = std::atan2(p.y(),p.x()) ;
658
659 if ( pPhi < -kAngTolerance*0.5 ) { pPhi += twopi ; } // 0<=pPhi<2pi
660 if ( fSPhi >= 0 )
661 {
662 if ( (std::abs(pPhi) < kAngTolerance*0.5)
663 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
664 {
665 pPhi += twopi ; // 0 <= pPhi < 2pi
666 }
667 if ( (pPhi >= fSPhi - kAngTolerance*0.5)
668 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
669 {
670 in = kSurface;
671 }
672 }
673 else // fSPhi < 0
674 {
675 if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
676 && (pPhi >= fSPhi + fDPhi + kAngTolerance*0.5) ) {;}
677 else
678 {
679 in = kSurface ;
680 }
681 }
682 }
683 }
684 }
685 return in ;
686}
687
688/////////////////////////////////////////////////////////////////////////////
689//
690// Return unit normal of surface closest to p
691// - note if point on z axis, ignore phi divided sides
692// - unsafe if point close to z axis a rmin=0 - no explicit checks
693
694G4ThreeVector G4Torus::SurfaceNormal( const G4ThreeVector& p ) const
695{
696 G4int noSurfaces = 0;
697 G4double rho2, rho, pt2, pt, pPhi;
698 G4double distRMin = kInfinity;
699 G4double distSPhi = kInfinity, distEPhi = kInfinity;
700 G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
701 G4ThreeVector nR, nPs, nPe;
702 G4ThreeVector norm, sumnorm(0.,0.,0.);
703
704 rho2 = p.x()*p.x() + p.y()*p.y();
705 rho = std::sqrt(rho2);
706 pt2 = std::fabs(rho2+p.z()*p.z() +fRtor*fRtor - 2*fRtor*rho);
707 pt = std::sqrt(pt2) ;
708
709 G4double distRMax = std::fabs(pt - fRmax);
710 if(fRmin) distRMin = std::fabs(pt - fRmin);
711
712 if( rho > delta )
713 {
714 nR = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
715 p.y()*(1-fRtor/rho)/pt,
716 p.z()/pt );
717 }
718
719 if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
720 {
721 if ( rho )
722 {
723 pPhi = std::atan2(p.y(),p.x());
724
725 if(pPhi < fSPhi-delta) { pPhi += twopi; }
726 else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
727
728 distSPhi = std::fabs( pPhi - fSPhi );
729 distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
730 }
731 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
732 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
733 }
734 if( distRMax <= delta )
735 {
736 noSurfaces ++;
737 sumnorm += nR;
738 }
739 if( fRmin && distRMin <= delta )
740 {
741 noSurfaces ++;
742 sumnorm -= nR;
743 }
744 if( fDPhi < twopi )
745 {
746 if (distSPhi <= dAngle)
747 {
748 noSurfaces ++;
749 sumnorm += nPs;
750 }
751 if (distEPhi <= dAngle)
752 {
753 noSurfaces ++;
754 sumnorm += nPe;
755 }
756 }
757 if ( noSurfaces == 0 )
758 {
759#ifdef G4CSGDEBUG
760 G4Exception("G4Torus::SurfaceNormal(p)", "Notification", JustWarning,
761 "Point p is not on surface !?" );
762#endif
763 norm = ApproxSurfaceNormal(p);
764 }
765 else if ( noSurfaces == 1 ) { norm = sumnorm; }
766 else { norm = sumnorm.unit(); }
767
768 return norm ;
769}
770
771//////////////////////////////////////////////////////////////////////////////
772//
773// Algorithm for SurfaceNormal() following the original specification
774// for points not on the surface
775
776G4ThreeVector G4Torus::ApproxSurfaceNormal( const G4ThreeVector& p ) const
777{
778 ENorm side ;
779 G4ThreeVector norm;
780 G4double rho2,rho,pt2,pt,phi;
781 G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
782
783 rho2 = p.x()*p.x() + p.y()*p.y();
784 rho = std::sqrt(rho2) ;
785 pt2 = std::fabs(rho2+p.z()*p.z() +fRtor*fRtor - 2*fRtor*rho) ;
786 pt = std::sqrt(pt2) ;
787
788 distRMax = std::fabs(pt - fRmax) ;
789
790 if(fRmin) // First minimum radius
791 {
792 distRMin = std::fabs(pt - fRmin) ;
793
794 if (distRMin < distRMax)
795 {
796 distMin = distRMin ;
797 side = kNRMin ;
798 }
799 else
800 {
801 distMin = distRMax ;
802 side = kNRMax ;
803 }
804 }
805 else
806 {
807 distMin = distRMax ;
808 side = kNRMax ;
809 }
810 if ( (fDPhi < twopi) && rho )
811 {
812 phi = std::atan2(p.y(),p.x()) ; // Protected against (0,0,z) (above rho!=0)
813
814 if (phi < 0) { phi += twopi ; }
815
816 if (fSPhi < 0 ) { distSPhi = std::fabs(phi-(fSPhi+twopi))*rho ; }
817 else { distSPhi = std::fabs(phi-fSPhi)*rho ; }
818
819 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
820
821 if (distSPhi < distEPhi) // Find new minimum
822 {
823 if (distSPhi<distMin) side = kNSPhi ;
824 }
825 else
826 {
827 if (distEPhi < distMin) { side = kNEPhi ; }
828 }
829 }
830 switch (side)
831 {
832 case kNRMin: // Inner radius
833 norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt,
834 -p.y()*(1-fRtor/rho)/pt,
835 -p.z()/pt ) ;
836 break ;
837 case kNRMax: // Outer radius
838 norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
839 p.y()*(1-fRtor/rho)/pt,
840 p.z()/pt ) ;
841 break;
842 case kNSPhi:
843 norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
844 break;
845 case kNEPhi:
846 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
847 break;
848 default:
849 DumpInfo();
850 G4Exception("G4Torus::ApproxSurfaceNormal()",
851 "Notification", JustWarning,
852 "Undefined side for valid surface normal to solid.");
853 break ;
854 }
855 return norm ;
856}
857
858///////////////////////////////////////////////////////////////////////
859//
860// Calculate distance to shape from outside, along normalised vector
861// - return kInfinity if no intersection, or intersection distance <= tolerance
862//
863// - Compute the intersection with the z planes
864// - if at valid r, phi, return
865//
866// -> If point is outer outer radius, compute intersection with rmax
867// - if at valid phi,z return
868//
869// -> Compute intersection with inner radius, taking largest +ve root
870// - if valid (phi), save intersction
871//
872// -> If phi segmented, compute intersections with phi half planes
873// - return smallest of valid phi intersections and
874// inner radius intersection
875//
876// NOTE:
877// - Precalculations for phi trigonometry are Done `just in time'
878// - `if valid' implies tolerant checking of intersection points
879
880G4double G4Torus::DistanceToIn( const G4ThreeVector& p,
881 const G4ThreeVector& v ) const
882{
883
884 G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value
885
886 G4double s[4] ;
887
888 // Precalculated trig for phi intersections - used by r,z intersections to
889 // check validity
890
891 G4bool seg; // true if segmented
892 G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT=0.,cosHDPhiIT=0.;
893 // half dphi + outer tolerance
894 G4double cPhi,sinCPhi=0.,cosCPhi=0.; // central phi
895
896 G4double tolORMin2,tolIRMin2; // `generous' radii squared
897 G4double tolORMax2,tolIRMax2 ;
898
899 G4double Dist,xi,yi,zi,rhoi2,it2; // Intersection point variables
900
901
902 G4double Comp;
903 G4double cosSPhi,sinSPhi; // Trig for phi start intersect
904 G4double ePhi,cosEPhi,sinEPhi; // for phi end intersect
905
906 // Set phi divided flag and precalcs
907 //
908 if ( fDPhi < twopi )
909 {
910 seg = true ;
911 hDPhi = 0.5*fDPhi ; // half delta phi
912 cPhi = fSPhi + hDPhi ;
913 hDPhiOT = hDPhi+0.5*kAngTolerance ; // outers tol' half delta phi
914 hDPhiIT = hDPhi - 0.5*kAngTolerance ;
915 sinCPhi = std::sin(cPhi) ;
916 cosCPhi = std::cos(cPhi) ;
917 cosHDPhiOT = std::cos(hDPhiOT) ;
918 cosHDPhiIT = std::cos(hDPhiIT) ;
919 }
920 else
921 {
922 seg = false ;
923 }
924
925 if (fRmin > kRadTolerance) // Calculate tolerant rmin and rmax
926 {
927 tolORMin2 = (fRmin - 0.5*kRadTolerance)*(fRmin - 0.5*kRadTolerance) ;
928 tolIRMin2 = (fRmin + 0.5*kRadTolerance)*(fRmin + 0.5*kRadTolerance) ;
929 }
930 else
931 {
932 tolORMin2 = 0 ;
933 tolIRMin2 = 0 ;
934 }
935 tolORMax2 = (fRmax + 0.5*kRadTolerance)*(fRmax + 0.5*kRadTolerance) ;
936 tolIRMax2 = (fRmax - kRadTolerance*0.5)*(fRmax - kRadTolerance*0.5) ;
937
938 // Intersection with Rmax (possible return) and Rmin (must also check phi)
939
940 G4double Rtor2 = fRtor*fRtor ;
941
942 snxt = SolveNumericJT(p,v,fRmax,true);
943 if (fRmin) // Possible Rmin intersection
944 {
945 s[0] = SolveNumericJT(p,v,fRmin,true);
946 if ( s[0] < snxt ) { snxt = s[0] ; }
947 }
948
949 //
950 // Phi segment intersection
951 //
952 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
953 //
954 // o NOTE: Large duplication of code between sphi & ephi checks
955 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
956 // intersection check <=0 -> >=0
957 // -> use some form of loop Construct ?
958
959 if (seg)
960 {
961 sinSPhi = std::sin(fSPhi) ; // First phi surface (`S'tarting phi)
962 cosSPhi = std::cos(fSPhi) ;
963 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; // Component in outwards
964 // normal direction
965 if (Comp < 0 )
966 {
967 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
968
969 if (Dist < kCarTolerance*0.5)
970 {
971 sphi = Dist/Comp ;
972 if (sphi < snxt)
973 {
974 if ( sphi < 0 ) { sphi = 0 ; }
975
976 xi = p.x() + sphi*v.x() ;
977 yi = p.y() + sphi*v.y() ;
978 zi = p.z() + sphi*v.z() ;
979 rhoi2 = xi*xi + yi*yi ;
980 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
981
982 if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
983 {
984 // r intersection is good - check intersecting
985 // with correct half-plane
986 //
987 if ((yi*cosCPhi-xi*sinCPhi)<=0) { snxt=sphi; }
988 }
989 }
990 }
991 }
992 ePhi=fSPhi+fDPhi; // Second phi surface (`E'nding phi)
993 sinEPhi=std::sin(ePhi);
994 cosEPhi=std::cos(ePhi);
995 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
996
997 if ( Comp < 0 ) // Component in outwards normal dirn
998 {
999 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1000
1001 if (Dist < kCarTolerance*0.5 )
1002 {
1003 sphi = Dist/Comp ;
1004 if (sphi < snxt )
1005 {
1006 if (sphi < 0 ) { sphi = 0 ; }
1007
1008 xi = p.x() + sphi*v.x() ;
1009 yi = p.y() + sphi*v.y() ;
1010 zi = p.z() + sphi*v.z() ;
1011 rhoi2 = xi*xi + yi*yi ;
1012 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
1013
1014 if (it2 >= tolORMin2 && it2 <= tolORMax2)
1015 {
1016 // z and r intersections good - check intersecting
1017 // with correct half-plane
1018 //
1019 if ((yi*cosCPhi-xi*sinCPhi)>=0) { snxt=sphi; }
1020 }
1021 }
1022 }
1023 }
1024 }
1025 if(snxt < 0.5*kCarTolerance) { snxt = 0.0 ; }
1026
1027 return snxt ;
1028}
1029
1030/////////////////////////////////////////////////////////////////////////////
1031//
1032// Calculate distance (<= actual) to closest surface of shape from outside
1033// - Calculate distance to z, radial planes
1034// - Only to phi planes if outside phi extent
1035// - Return 0 if point inside
1036
1037G4double G4Torus::DistanceToIn( const G4ThreeVector& p ) const
1038{
1039 G4double safe=0.0, safe1, safe2 ;
1040 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1041 G4double rho2, rho, pt2, pt ;
1042
1043 rho2 = p.x()*p.x() + p.y()*p.y() ;
1044 rho = std::sqrt(rho2) ;
1045 pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
1046 pt = std::sqrt(pt2) ;
1047
1048 safe1 = fRmin - pt ;
1049 safe2 = pt - fRmax ;
1050
1051 if (safe1 > safe2) { safe = safe1; }
1052 else { safe = safe2; }
1053
1054 if ( fDPhi < twopi && rho )
1055 {
1056 phiC = fSPhi + fDPhi*0.5 ;
1057 cosPhiC = std::cos(phiC) ;
1058 sinPhiC = std::sin(phiC) ;
1059 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1060
1061 if (cosPsi < std::cos(fDPhi*0.5) ) // Psi=angle from central phi to point
1062 { // Point lies outside phi range
1063 if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1064 {
1065 safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1066 }
1067 else
1068 {
1069 ePhi = fSPhi + fDPhi ;
1070 safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1071 }
1072 if (safePhi > safe) { safe = safePhi ; }
1073 }
1074 }
1075 if (safe < 0 ) { safe = 0 ; }
1076 return safe;
1077}
1078
1079///////////////////////////////////////////////////////////////////////////
1080//
1081// Calculate distance to surface of shape from `inside', allowing for tolerance
1082// - Only Calc rmax intersection if no valid rmin intersection
1083//
1084
1085G4double G4Torus::DistanceToOut( const G4ThreeVector& p,
1086 const G4ThreeVector& v,
1087 const G4bool calcNorm,
1088 G4bool *validNorm,
1089 G4ThreeVector *n ) const
1090{
1091 ESide side = kNull, sidephi = kNull ;
1092 G4double snxt = kInfinity, sphi, s[4] ;
1093
1094 // Vars for phi intersection
1095 //
1096 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
1097 G4double cPhi, sinCPhi, cosCPhi ;
1098 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
1099
1100 // Radial Intersections Defenitions & General Precals
1101
1102 //////////////////////// new calculation //////////////////////
1103
1104#if 1
1105
1106 // This is the version with the calculation of CalcNorm = true
1107 // To be done: Check the precision of this calculation.
1108 // If you want return always validNorm = false, then take the version below
1109
1110 G4double Rtor2 = fRtor*fRtor ;
1111 G4double rho2 = p.x()*p.x()+p.y()*p.y();
1112 G4double rho = std::sqrt(rho2) ;
1113
1114
1115 G4double pt2 = std::fabs(rho2 + p.z()*p.z() + Rtor2 - 2*fRtor*rho) ;
1116 G4double pt = std::sqrt(pt2) ;
1117
1118 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
1119
1120 G4double tolRMax = fRmax - kRadTolerance*0.5 ;
1121
1122 G4double vDotNmax = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
1123 G4double pDotxyNmax = (1 - fRtor/rho) ;
1124
1125 if( (pt2 > tolRMax*tolRMax) && (vDotNmax >= 0) )
1126 {
1127 // On tolerant boundary & heading outwards (or perpendicular to) outer
1128 // radial surface -> leaving immediately with *n for really convex part
1129 // only
1130
1131 if ( calcNorm && (pDotxyNmax >= -kRadTolerance) )
1132 {
1133 *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
1134 p.y()*(1 - fRtor/rho)/pt,
1135 p.z()/pt ) ;
1136 *validNorm = true ;
1137 }
1138 return snxt = 0 ; // Leaving by Rmax immediately
1139 }
1140
1141 snxt = SolveNumericJT(p,v,fRmax,false);
1142 side = kRMax ;
1143
1144 // rmin
1145
1146 if ( fRmin )
1147 {
1148 G4double tolRMin = fRmin + kRadTolerance*0.5 ;
1149
1150 if ( (pt2 < tolRMin*tolRMin) && (vDotNmax < 0) )
1151 {
1152 if (calcNorm) { *validNorm = false ; } // Concave surface of the torus
1153 return snxt = 0 ; // Leaving by Rmin immediately
1154 }
1155
1156 s[0] = SolveNumericJT(p,v,fRmin,false);
1157 if ( s[0] < snxt )
1158 {
1159 snxt = s[0] ;
1160 side = kRMin ;
1161 }
1162 }
1163
1164#else
1165
1166 // this is the "conservative" version which return always validnorm = false
1167 // NOTE: using this version the unit test testG4Torus will break
1168
1169 snxt = SolveNumericJT(p,v,fRmax,false);
1170 side = kRMax ;
1171
1172 if ( fRmin )
1173 {
1174 s[0] = SolveNumericJT(p,v,fRmin,false);
1175 if ( s[0] < snxt )
1176 {
1177 snxt = s[0] ;
1178 side = kRMin ;
1179 }
1180 }
1181
1182 if ( calcNorm && (snxt == 0.0) )
1183 {
1184 *validNorm = false ; // Leaving solid, but possible re-intersection
1185 return snxt ;
1186 }
1187
1188#endif
1189
1190 if (fDPhi < twopi) // Phi Intersections
1191 {
1192 sinSPhi = std::sin(fSPhi) ;
1193 cosSPhi = std::cos(fSPhi) ;
1194 ePhi = fSPhi + fDPhi ;
1195 sinEPhi = std::sin(ePhi) ;
1196 cosEPhi = std::cos(ePhi) ;
1197 cPhi = fSPhi + fDPhi*0.5 ;
1198 sinCPhi = std::sin(cPhi) ;
1199 cosCPhi = std::cos(cPhi) ;
1200
1201 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1202 {
1203 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside
1204 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1205
1206 // Comp -ve when in direction of outwards normal
1207 //
1208 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1209 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1210 sidephi = kNull ;
1211
1212 if ( (pDistS <= 0) && (pDistE <= 0) )
1213 {
1214 // Inside both phi *full* planes
1215
1216 if (compS<0)
1217 {
1218 sphi=pDistS/compS;
1219 xi=p.x()+sphi*v.x();
1220 yi=p.y()+sphi*v.y();
1221
1222 // Check intersecting with correct half-plane
1223 // (if not -> no intersect)
1224 //
1225 if ((yi*cosCPhi-xi*sinCPhi)>=0)
1226 {
1227 sphi=kInfinity;
1228 }
1229 else
1230 {
1231 sidephi=kSPhi;
1232 if (pDistS>-kCarTolerance*0.5) { sphi=0; } // Leave by sphi
1233 // immediately
1234 }
1235 }
1236 else
1237 {
1238 sphi=kInfinity;
1239 }
1240
1241 if (compE<0)
1242 {
1243 sphi2=pDistE/compE;
1244
1245 // Only check further if < starting phi intersection
1246 //
1247 if (sphi2<sphi)
1248 {
1249 xi=p.x()+sphi2*v.x();
1250 yi=p.y()+sphi2*v.y();
1251
1252 // Check intersecting with correct half-plane
1253 //
1254 if ((yi*cosCPhi-xi*sinCPhi)>=0)
1255 {
1256 // Leaving via ending phi
1257 //
1258 sidephi=kEPhi;
1259 if (pDistE<=-kCarTolerance*0.5)
1260 {
1261 sphi=sphi2;
1262 }
1263 else
1264 {
1265 sphi=0;
1266 }
1267 }
1268 }
1269 }
1270 }
1271 else if ( (pDistS>=0) && (pDistE>=0) )
1272 {
1273 // Outside both *full* phi planes
1274
1275 if (pDistS <= pDistE)
1276 {
1277 sidephi = kSPhi ;
1278 }
1279 else
1280 {
1281 sidephi = kEPhi ;
1282 }
1283 if (fDPhi>pi)
1284 {
1285 if ( (compS<0) && (compE<0) ) { sphi=0; }
1286 else { sphi=kInfinity; }
1287 }
1288 else
1289 {
1290 // if towards both >=0 then once inside (after error)
1291 // will remain inside
1292 //
1293 if ( (compS>=0) && (compE>=0) )
1294 {
1295 sphi=kInfinity;
1296 }
1297 else
1298 {
1299 sphi=0;
1300 }
1301 }
1302 }
1303 else if ( (pDistS>0) && (pDistE<0) )
1304 {
1305 // Outside full starting plane, inside full ending plane
1306
1307 if (fDPhi>pi)
1308 {
1309 if (compE<0)
1310 {
1311 sphi=pDistE/compE;
1312 xi=p.x()+sphi*v.x();
1313 yi=p.y()+sphi*v.y();
1314
1315 // Check intersection in correct half-plane
1316 // (if not -> not leaving phi extent)
1317 //
1318 if ((yi*cosCPhi-xi*sinCPhi)<=0)
1319 {
1320 sphi=kInfinity;
1321 }
1322 else
1323 {
1324 // Leaving via Ending phi
1325 //
1326 sidephi = kEPhi ;
1327 if (pDistE>-kCarTolerance*0.5) { sphi=0; }
1328 }
1329 }
1330 else
1331 {
1332 sphi=kInfinity;
1333 }
1334 }
1335 else
1336 {
1337 if (compS>=0)
1338 {
1339 if (compE<0)
1340 {
1341 sphi=pDistE/compE;
1342 xi=p.x()+sphi*v.x();
1343 yi=p.y()+sphi*v.y();
1344
1345 // Check intersection in correct half-plane
1346 // (if not -> remain in extent)
1347 //
1348 if ((yi*cosCPhi-xi*sinCPhi)<=0)
1349 {
1350 sphi=kInfinity;
1351 }
1352 else
1353 {
1354 // otherwise leaving via Ending phi
1355 //
1356 sidephi=kEPhi;
1357 }
1358 }
1359 else { sphi=kInfinity; }
1360 }
1361 else
1362 {
1363 // leaving immediately by starting phi
1364 //
1365 sidephi=kSPhi;
1366 sphi=0;
1367 }
1368 }
1369 }
1370 else
1371 {
1372 // Must be pDistS<0&&pDistE>0
1373 // Inside full starting plane, outside full ending plane
1374
1375 if (fDPhi>pi)
1376 {
1377 if (compS<0)
1378 {
1379 sphi=pDistS/compS;
1380 xi=p.x()+sphi*v.x();
1381 yi=p.y()+sphi*v.y();
1382
1383 // Check intersection in correct half-plane
1384 // (if not -> not leaving phi extent)
1385 //
1386 if ((yi*cosCPhi-xi*sinCPhi)>=0)
1387 {
1388 sphi=kInfinity;
1389 }
1390 else
1391 {
1392 // Leaving via Starting phi
1393 //
1394 sidephi = kSPhi ;
1395 if (pDistS>-kCarTolerance*0.5) { sphi=0; }
1396 }
1397 }
1398 else
1399 {
1400 sphi=kInfinity;
1401 }
1402 }
1403 else
1404 {
1405 if (compE>=0)
1406 {
1407 if (compS<0)
1408 {
1409 sphi=pDistS/compS;
1410 xi=p.x()+sphi*v.x();
1411 yi=p.y()+sphi*v.y();
1412
1413 // Check intersection in correct half-plane
1414 // (if not -> remain in extent)
1415 //
1416 if ((yi*cosCPhi-xi*sinCPhi)>=0)
1417 {
1418 sphi=kInfinity;
1419 }
1420 else
1421 {
1422 // otherwise leaving via Starting phi
1423 //
1424 sidephi=kSPhi;
1425 }
1426 }
1427 else { sphi=kInfinity; }
1428 }
1429 else
1430 {
1431 // leaving immediately by ending
1432 //
1433 sidephi=kEPhi;
1434 sphi=0;
1435 }
1436 }
1437 }
1438 }
1439 else
1440 {
1441 // On z axis + travel not || to z axis -> if phi of vector direction
1442 // within phi of shape, Step limited by rmax, else Step =0
1443
1444 vphi=std::atan2(v.y(),v.x());
1445 if ( (fSPhi<vphi) && (vphi<fSPhi+fDPhi) )
1446 {
1447 sphi=kInfinity;
1448 }
1449 else
1450 {
1451 sidephi = kSPhi ; // arbitrary
1452 sphi=0;
1453 }
1454 }
1455
1456 // Order intersections
1457
1458 if (sphi<snxt)
1459 {
1460 snxt=sphi;
1461 side=sidephi;
1462 }
1463 }
1464 G4double rhoi2,rhoi,it2,it,iDotxyNmax ;
1465
1466 // Note: by numerical computation we know where the ray hits the torus
1467 // So I propose to return the side where the ray hits
1468
1469 if (calcNorm)
1470 {
1471 switch(side)
1472 {
1473 case kRMax: // n is unit vector
1474 xi = p.x() + snxt*v.x() ;
1475 yi =p.y() + snxt*v.y() ;
1476 zi = p.z() + snxt*v.z() ;
1477 rhoi2 = xi*xi + yi*yi ;
1478 rhoi = std::sqrt(rhoi2) ;
1479 it2 = std::fabs(rhoi2 + zi*zi + fRtor*fRtor - 2*fRtor*rhoi) ;
1480 it = std::sqrt(it2) ;
1481 iDotxyNmax = (1-fRtor/rhoi) ;
1482 if(iDotxyNmax >= -kRadTolerance) // really convex part of Rmax
1483 {
1484 *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
1485 yi*(1-fRtor/rhoi)/it,
1486 zi/it ) ;
1487 *validNorm = true ;
1488 }
1489 else
1490 {
1491 *validNorm = false ; // concave-convex part of Rmax
1492 }
1493 break ;
1494
1495 case kRMin:
1496 *validNorm = false ; // Rmin is concave or concave-convex
1497 break;
1498
1499 case kSPhi:
1500 if (fDPhi <= pi )
1501 {
1502 *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
1503 *validNorm=true;
1504 }
1505 else
1506 {
1507 *validNorm = false ;
1508 }
1509 break ;
1510
1511 case kEPhi:
1512 if (fDPhi <= pi)
1513 {
1514 *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
1515 *validNorm=true;
1516 }
1517 else
1518 {
1519 *validNorm = false ;
1520 }
1521 break;
1522
1523 default:
1524
1525 // It seems we go here from time to time ...
1526
1527 G4cout.precision(16);
1528 G4cout << G4endl;
1529 DumpInfo();
1530 G4cout << "Position:" << G4endl << G4endl;
1531 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
1532 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
1533 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
1534 G4cout << "Direction:" << G4endl << G4endl;
1535 G4cout << "v.x() = " << v.x() << G4endl;
1536 G4cout << "v.y() = " << v.y() << G4endl;
1537 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
1538 G4cout << "Proposed distance :" << G4endl << G4endl;
1539 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl;
1540 G4Exception("G4Torus::DistanceToOut(p,v,..)",
1541 "Notification",JustWarning,
1542 "Undefined side for valid surface normal to solid.");
1543 break;
1544 }
1545 }
1546
1547 return snxt;
1548}
1549
1550/////////////////////////////////////////////////////////////////////////
1551//
1552// Calculate distance (<=actual) to closest surface of shape from inside
1553
1554G4double G4Torus::DistanceToOut( const G4ThreeVector& p ) const
1555{
1556 G4double safe=0.0,safeR1,safeR2;
1557 G4double rho2,rho,pt2,pt ;
1558 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1559 rho2 = p.x()*p.x() + p.y()*p.y() ;
1560 rho = std::sqrt(rho2) ;
1561 pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
1562 pt = std::sqrt(pt2) ;
1563
1564#ifdef G4CSGDEBUG
1565 if( Inside(p) == kOutside )
1566 {
1567 G4cout.precision(16) ;
1568 G4cout << G4endl ;
1569 DumpInfo();
1570 G4cout << "Position:" << G4endl << G4endl ;
1571 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1572 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1573 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1574 G4Exception("G4Torus::DistanceToOut(p)", "Notification",
1575 JustWarning, "Point p is outside !?" );
1576 }
1577#endif
1578
1579 if (fRmin)
1580 {
1581 safeR1 = pt - fRmin ;
1582 safeR2 = fRmax - pt ;
1583
1584 if (safeR1 < safeR2) { safe = safeR1 ; }
1585 else { safe = safeR2 ; }
1586 }
1587 else
1588 {
1589 safe = fRmax - pt ;
1590 }
1591
1592 // Check if phi divided, Calc distances closest phi plane
1593 //
1594 if (fDPhi<twopi) // Above/below central phi of Torus?
1595 {
1596 phiC = fSPhi + fDPhi*0.5 ;
1597 cosPhiC = std::cos(phiC) ;
1598 sinPhiC = std::sin(phiC) ;
1599
1600 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1601 {
1602 safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1603 }
1604 else
1605 {
1606 ePhi = fSPhi + fDPhi ;
1607 safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1608 }
1609 if (safePhi < safe) { safe = safePhi ; }
1610 }
1611 if (safe < 0) { safe = 0 ; }
1612 return safe ;
1613}
1614
1615/////////////////////////////////////////////////////////////////////////////
1616//
1617// Create a List containing the transformed vertices
1618// Ordering [0-3] -fRtor cross section
1619// [4-7] +fRtor cross section such that [0] is below [4],
1620// [1] below [5] etc.
1621// Note:
1622// Caller has deletion resposibility
1623// Potential improvement: For last slice, use actual ending angle
1624// to avoid rounding error problems.
1625
1626G4ThreeVectorList*
1627G4Torus::CreateRotatedVertices( const G4AffineTransform& pTransform,
1628 G4int& noPolygonVertices ) const
1629{
1630 G4ThreeVectorList *vertices;
1631 G4ThreeVector vertex0,vertex1,vertex2,vertex3;
1632 G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle;
1633 G4double rMaxX,rMaxY,rMinX,rMinY;
1634 G4int crossSection,noCrossSections;
1635
1636 // Compute no of cross-sections necessary to mesh tube
1637 //
1638 noCrossSections = G4int (fDPhi/kMeshAngleDefault) + 1 ;
1639
1640 if (noCrossSections < kMinMeshSections)
1641 {
1642 noCrossSections = kMinMeshSections ;
1643 }
1644 else if (noCrossSections>kMaxMeshSections)
1645 {
1646 noCrossSections=kMaxMeshSections;
1647 }
1648 meshAngle = fDPhi/(noCrossSections - 1) ;
1649 meshRMax = (fRtor + fRmax)/std::cos(meshAngle*0.5) ;
1650
1651 // If complete in phi, set start angle such that mesh will be at fRmax
1652 // on the x axis. Will give better extent calculations when not rotated
1653
1654 if ( (fDPhi == pi*2.0) && (fSPhi == 0) )
1655 {
1656 sAngle = -meshAngle*0.5 ;
1657 }
1658 else
1659 {
1660 sAngle = fSPhi ;
1661 }
1662 vertices = new G4ThreeVectorList();
1663 vertices->reserve(noCrossSections*4) ;
1664
1665 if (vertices)
1666 {
1667 for (crossSection=0;crossSection<noCrossSections;crossSection++)
1668 {
1669 // Compute coordinates of cross section at section crossSection
1670
1671 crossAngle=sAngle+crossSection*meshAngle;
1672 cosCrossAngle=std::cos(crossAngle);
1673 sinCrossAngle=std::sin(crossAngle);
1674
1675 rMaxX=meshRMax*cosCrossAngle;
1676 rMaxY=meshRMax*sinCrossAngle;
1677 rMinX=(fRtor-fRmax)*cosCrossAngle;
1678 rMinY=(fRtor-fRmax)*sinCrossAngle;
1679 vertex0=G4ThreeVector(rMinX,rMinY,-fRmax);
1680 vertex1=G4ThreeVector(rMaxX,rMaxY,-fRmax);
1681 vertex2=G4ThreeVector(rMaxX,rMaxY,+fRmax);
1682 vertex3=G4ThreeVector(rMinX,rMinY,+fRmax);
1683
1684 vertices->push_back(pTransform.TransformPoint(vertex0));
1685 vertices->push_back(pTransform.TransformPoint(vertex1));
1686 vertices->push_back(pTransform.TransformPoint(vertex2));
1687 vertices->push_back(pTransform.TransformPoint(vertex3));
1688 }
1689 noPolygonVertices = 4 ;
1690 }
1691 else
1692 {
1693 DumpInfo();
1694 G4Exception("G4Torus::CreateRotatedVertices()",
1695 "FatalError", FatalException,
1696 "Error in allocation of vertices. Out of memory !");
1697 }
1698 return vertices;
1699}
1700
1701//////////////////////////////////////////////////////////////////////////
1702//
1703// Stream object contents to an output stream
1704
1705G4GeometryType G4Torus::GetEntityType() const
1706{
1707 return G4String("G4Torus");
1708}
1709
1710//////////////////////////////////////////////////////////////////////////
1711//
1712// Stream object contents to an output stream
1713
1714std::ostream& G4Torus::StreamInfo( std::ostream& os ) const
1715{
1716 os << "-----------------------------------------------------------\n"
1717 << " *** Dump for solid - " << GetName() << " ***\n"
1718 << " ===================================================\n"
1719 << " Solid type: G4Torus\n"
1720 << " Parameters: \n"
1721 << " inner radius: " << fRmin/mm << " mm \n"
1722 << " outer radius: " << fRmax/mm << " mm \n"
1723 << " swept radius: " << fRtor/mm << " mm \n"
1724 << " starting phi: " << fSPhi/degree << " degrees \n"
1725 << " delta phi : " << fDPhi/degree << " degrees \n"
1726 << "-----------------------------------------------------------\n";
1727
1728 return os;
1729}
1730
1731////////////////////////////////////////////////////////////////////////////
1732//
1733// GetPointOnSurface
1734
1735G4ThreeVector G4Torus::GetPointOnSurface() const
1736{
1737 G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
1738
1739 phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
1740 theta = RandFlat::shoot(0.,2.*pi);
1741
1742 cosu = std::cos(phi); sinu = std::sin(phi);
1743 cosv = std::cos(theta); sinv = std::sin(theta);
1744
1745 // compute the areas
1746
1747 aOut = (fDPhi)*2.*pi*fRtor*fRmax;
1748 aIn = (fDPhi)*2.*pi*fRtor*fRmin;
1749 aSide = pi*(fRmax*fRmax-fRmin*fRmin);
1750
1751 if(fSPhi == 0 && fDPhi == twopi){ aSide = 0; }
1752 chose = RandFlat::shoot(0.,aOut + aIn + 2.*aSide);
1753
1754 if(chose < aOut)
1755 {
1756 return G4ThreeVector ((fRtor+fRmax*cosv)*cosu,
1757 (fRtor+fRmax*cosv)*sinu, fRmax*sinv);
1758 }
1759 else if( (chose >= aOut) && (chose < aOut + aIn) )
1760 {
1761 return G4ThreeVector ((fRtor+fRmin*cosv)*cosu,
1762 (fRtor+fRmin*cosv)*sinu, fRmin*sinv);
1763 }
1764 else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
1765 {
1766 rRand = RandFlat::shoot(fRmin,fRmax);
1767 return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi),
1768 (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv);
1769 }
1770 else
1771 {
1772 rRand = RandFlat::shoot(fRmin,fRmax);
1773 return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
1774 (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi),
1775 rRand*sinv);
1776 }
1777}
1778
1779///////////////////////////////////////////////////////////////////////
1780//
1781// Visualisation Functions
1782
1783void G4Torus::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1784{
1785 scene.AddSolid (*this);
1786}
1787
1788G4Polyhedron* G4Torus::CreatePolyhedron () const
1789{
1790 return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi);
1791}
1792
1793G4NURBS* G4Torus::CreateNURBS () const
1794{
1795 G4NURBS* pNURBS;
1796 if (fRmin != 0)
1797 {
1798 if (fDPhi >= 2.0 * pi)
1799 {
1800 pNURBS = new G4NURBStube(fRmin, fRmax, fRtor);
1801 }
1802 else
1803 {
1804 pNURBS = new G4NURBStubesector(fRmin, fRmax, fRtor, fSPhi, fSPhi + fDPhi);
1805 }
1806 }
1807 else
1808 {
1809 if (fDPhi >= 2.0 * pi)
1810 {
1811 pNURBS = new G4NURBScylinder (fRmax, fRtor);
1812 }
1813 else
1814 {
1815 const G4double epsilon = 1.e-4; // Cylinder sector not yet available!
1816 pNURBS = new G4NURBStubesector (epsilon, fRmax, fRtor,
1817 fSPhi, fSPhi + fDPhi);
1818 }
1819 }
1820 return pNURBS;
1821}
Note: See TracBrowser for help on using the repository browser.