source: trunk/source/geometry/solids/CSG/src/G4Tubs.cc@ 1090

Last change on this file since 1090 was 1058, checked in by garnier, 17 years ago

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