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

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

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