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

Last change on this file since 1051 was 850, checked in by garnier, 17 years ago

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