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

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

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