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

Last change on this file since 833 was 831, checked in by garnier, 17 years ago

import all except CVS
File size: 55.1 KB
RevLine 
[831]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4Tubs.cc,v 1.66 2007/11/23 09:07:43 tnikitin Exp $
28// GEANT4 tag $Name: $
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 ;
985 if (cosPsi >= cosHDPhiIT) { return 0.0; }
986 }
987 else
988 {
989 return 0.0 ;
990 }
991 }
992 }
993 if ( fRMin ) // Try inner cylinder intersection
994 {
995 c = (t3 - fRMin*fRMin)/t1 ;
996 d = b*b - c ;
997 if ( d >= 0.0 ) // If real root
998 {
999 // Always want 2nd root - we are outside and know rmax Hit was bad
1000 // - If on surface of rmin also need farthest root
1001
1002 s = -b + std::sqrt(d) ;
1003 if (s >= -0.5*kCarTolerance) // check forwards
1004 {
1005 // Check z intersection
1006 //
1007 if(s < 0.0) { s = 0.0; }
1008 zi = p.z() + s*v.z() ;
1009 if (std::fabs(zi) <= tolODz)
1010 {
1011 // Z ok. Check phi
1012 //
1013 if ( !seg )
1014 {
1015 return s ;
1016 }
1017 else
1018 {
1019 xi = p.x() + s*v.x() ;
1020 yi = p.y() + s*v.y() ;
1021 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1022 if (cosPsi >= cosHDPhiIT)
1023 {
1024 // Good inner radius isect
1025 // - but earlier phi isect still possible
1026
1027 snxt = s ;
1028 }
1029 }
1030 } // end if std::fabs(zi)
1031 } // end if (s>=0)
1032 } // end if (d>=0)
1033 } // end if (fRMin)
1034 }
1035
1036 // Phi segment intersection
1037 //
1038 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1039 //
1040 // o NOTE: Large duplication of code between sphi & ephi checks
1041 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1042 // intersection check <=0 -> >=0
1043 // -> use some form of loop Construct ?
1044 //
1045 if ( seg )
1046 {
1047 // First phi surface (Starting phi)
1048
1049 sinSPhi = std::sin(fSPhi) ;
1050 cosSPhi = std::cos(fSPhi) ;
1051 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1052
1053 if ( Comp < 0 ) // Component in outwards normal dirn
1054 {
1055 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1056
1057 if ( Dist < kCarTolerance*0.5 )
1058 {
1059 s = Dist/Comp ;
1060
1061 if (s < snxt)
1062 {
1063 if ( s < 0 ) { s = 0.0; }
1064 zi = p.z() + s*v.z() ;
1065 if ( std::fabs(zi) <= tolODz )
1066 {
1067 xi = p.x() + s*v.x() ;
1068 yi = p.y() + s*v.y() ;
1069 rho2 = xi*xi + yi*yi ;
1070
1071 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1072 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1073 && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 )
1074 && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) )
1075 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1076 && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1077 && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) )
1078 {
1079 // z and r intersections good
1080 // - check intersecting with correct half-plane
1081 //
1082 if ((yi*cosCPhi-xi*sinCPhi) <= 0) snxt = s ;
1083 }
1084 }
1085 }
1086 }
1087 }
1088
1089 // Second phi surface (`E'nding phi)
1090
1091 ePhi = fSPhi + fDPhi ;
1092 sinEPhi = std::sin(ePhi) ;
1093 cosEPhi = std::cos(ePhi) ;
1094 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1095
1096 if (Comp < 0 ) // Component in outwards normal dirn
1097 {
1098 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1099
1100 if ( Dist < kCarTolerance*0.5 )
1101 {
1102 s = Dist/Comp ;
1103
1104 if (s < snxt)
1105 {
1106 if ( s < 0 ) { s = 0; }
1107 zi = p.z() + s*v.z() ;
1108 if ( std::fabs(zi) <= tolODz )
1109 {
1110 xi = p.x() + s*v.x() ;
1111 yi = p.y() + s*v.y() ;
1112 rho2 = xi*xi + yi*yi ;
1113 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1114 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1115 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1116 && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1117 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1118 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1119 && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1120 {
1121 // z and r intersections good
1122 // - check intersecting with correct half-plane
1123 //
1124 if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = s; }
1125 }
1126 }
1127 }
1128 }
1129 } // Comp < 0
1130 } // seg != 0
1131 if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
1132 return snxt ;
1133}
1134
1135//////////////////////////////////////////////////////////////////
1136//
1137// Calculate distance to shape from outside, along normalised vector
1138// - return kInfinity if no intersection, or intersection distance <= tolerance
1139//
1140// - Compute the intersection with the z planes
1141// - if at valid r, phi, return
1142//
1143// -> If point is outer outer radius, compute intersection with rmax
1144// - if at valid phi,z return
1145//
1146// -> Compute intersection with inner radius, taking largest +ve root
1147// - if valid (in z,phi), save intersction
1148//
1149// -> If phi segmented, compute intersections with phi half planes
1150// - return smallest of valid phi intersections and
1151// inner radius intersection
1152//
1153// NOTE:
1154// - Precalculations for phi trigonometry are Done `just in time'
1155// - `if valid' implies tolerant checking of intersection points
1156// Calculate distance (<= actual) to closest surface of shape from outside
1157// - Calculate distance to z, radial planes
1158// - Only to phi planes if outside phi extent
1159// - Return 0 if point inside
1160
1161G4double G4Tubs::DistanceToIn( const G4ThreeVector& p ) const
1162{
1163 G4double safe=0.0, rho, safe1, safe2, safe3 ;
1164 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1165
1166 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1167 safe1 = fRMin - rho ;
1168 safe2 = rho - fRMax ;
1169 safe3 = std::fabs(p.z()) - fDz ;
1170
1171 if ( safe1 > safe2 ) { safe = safe1; }
1172 else { safe = safe2; }
1173 if ( safe3 > safe ) { safe = safe3; }
1174
1175 if (fDPhi < twopi && rho)
1176 {
1177 phiC = fSPhi + fDPhi*0.5 ;
1178 cosPhiC = std::cos(phiC) ;
1179 sinPhiC = std::sin(phiC) ;
1180
1181 // Psi=angle from central phi to point
1182 //
1183 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1184
1185 if ( cosPsi < std::cos(fDPhi*0.5) )
1186 {
1187 // Point lies outside phi range
1188
1189 if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1190 {
1191 safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1192 }
1193 else
1194 {
1195 ePhi = fSPhi + fDPhi ;
1196 safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1197 }
1198 if ( safePhi > safe ) { safe = safePhi; }
1199 }
1200 }
1201 if ( safe < 0 ) { safe = 0; }
1202 return safe ;
1203}
1204
1205//////////////////////////////////////////////////////////////////////////////
1206//
1207// Calculate distance to surface of shape from `inside', allowing for tolerance
1208// - Only Calc rmax intersection if no valid rmin intersection
1209
1210G4double G4Tubs::DistanceToOut( const G4ThreeVector& p,
1211 const G4ThreeVector& v,
1212 const G4bool calcNorm,
1213 G4bool *validNorm,
1214 G4ThreeVector *n ) const
1215{
1216 ESide side = kNull , sider = kNull, sidephi = kNull ;
1217 G4double snxt, sr = kInfinity, sphi = kInfinity, pdist ;
1218 G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1219
1220 // Vars for phi intersection:
1221
1222 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;
1223 G4double cPhi, sinCPhi, cosCPhi ;
1224 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1225
1226 // Z plane intersection
1227
1228 if (v.z() > 0 )
1229 {
1230 pdist = fDz - p.z() ;
1231 if ( pdist > kCarTolerance*0.5 )
1232 {
1233 snxt = pdist/v.z() ;
1234 side = kPZ ;
1235 }
1236 else
1237 {
1238 if (calcNorm)
1239 {
1240 *n = G4ThreeVector(0,0,1) ;
1241 *validNorm = true ;
1242 }
1243 return snxt = 0 ;
1244 }
1245 }
1246 else if ( v.z() < 0 )
1247 {
1248 pdist = fDz + p.z() ;
1249
1250 if ( pdist > kCarTolerance*0.5 )
1251 {
1252 snxt = -pdist/v.z() ;
1253 side = kMZ ;
1254 }
1255 else
1256 {
1257 if (calcNorm)
1258 {
1259 *n = G4ThreeVector(0,0,-1) ;
1260 *validNorm = true ;
1261 }
1262 return snxt = 0.0 ;
1263 }
1264 }
1265 else
1266 {
1267 snxt = kInfinity ; // Travel perpendicular to z axis
1268 side = kNull;
1269 }
1270
1271 // Radial Intersections
1272 //
1273 // Find intersction with cylinders at rmax/rmin
1274 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1275 //
1276 // Intersects with x^2+y^2=R^2
1277 //
1278 // 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
1279 //
1280 // t1 t2 t3
1281
1282 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1283 t2 = p.x()*v.x() + p.y()*v.y() ;
1284 t3 = p.x()*p.x() + p.y()*p.y() ;
1285
1286 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; }
1287 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz
1288
1289 if ( t1 > 0 ) // Check not parallel
1290 {
1291 // Calculate sr, r exit distance
1292
1293 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1294 {
1295 // Delta r not negative => leaving via rmax
1296
1297 deltaR = t3 - fRMax*fRMax ;
1298
1299 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1300 // - avoid sqrt for efficiency
1301
1302 if ( deltaR < -kRadTolerance*fRMax )
1303 {
1304 b = t2/t1 ;
1305 c = deltaR/t1 ;
1306 d2= b*b-c;
1307 if(d2>=0.){sr = -b + std::sqrt(d2);}
1308 else{sr=0.;};
1309 sider = kRMax ;
1310 }
1311 else
1312 {
1313 // On tolerant boundary & heading outwards (or perpendicular to)
1314 // outer radial surface -> leaving immediately
1315
1316 if ( calcNorm )
1317 {
1318 // if ( p.x() || p.y() )
1319 // {
1320 // *n=G4ThreeVector(p.x(),p.y(),0);
1321 // }
1322 // else
1323 // {
1324 // *n=v;
1325 // }
1326 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1327 *validNorm = true ;
1328 }
1329 return snxt = 0 ; // Leaving by rmax immediately
1330 }
1331 }
1332 else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection
1333 {
1334 roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1335
1336 if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1337 {
1338 deltaR = t3 - fRMin*fRMin ;
1339 b = t2/t1 ;
1340 c = deltaR/t1 ;
1341 d2 = b*b - c ;
1342
1343 if ( d2 >= 0 ) // Leaving via rmin
1344 {
1345 // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1346 // - avoid sqrt for efficiency
1347
1348 if (deltaR > kRadTolerance*fRMin)
1349 {
1350 sr = -b-std::sqrt(d2) ;
1351 sider = kRMin ;
1352 }
1353 else
1354 {
1355 if ( calcNorm ) { *validNorm = false; } // Concave side
1356 return snxt = 0.0;
1357 }
1358 }
1359 else // No rmin intersect -> must be rmax intersect
1360 {
1361 deltaR = t3 - fRMax*fRMax ;
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 perpendiculair 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 else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1382 // No rmin intersect -> must be rmax intersect
1383 {
1384 deltaR = t3 - fRMax*fRMax ;
1385 b = t2/t1 ;
1386 c = deltaR/t1;
1387 d2 = b*b-c;
1388 if(d2>=0.)
1389 {
1390 sr = -b + std::sqrt(d2) ;
1391 sider = kRMax ;
1392 }
1393 else // Case: On the border+t2<kRadTolerance
1394 // (v is perpendiculair to the surface)
1395 {
1396 if (calcNorm)
1397 {
1398 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1399 *validNorm = true ;
1400 }
1401 return snxt = 0.0;
1402 }
1403 }
1404 }
1405
1406 // Phi Intersection
1407
1408 if ( fDPhi < twopi )
1409 {
1410 sinSPhi = std::sin(fSPhi) ;
1411 cosSPhi = std::cos(fSPhi) ;
1412 ePhi = fSPhi + fDPhi ;
1413 sinEPhi = std::sin(ePhi) ;
1414 cosEPhi = std::cos(ePhi) ;
1415 cPhi = fSPhi + fDPhi*0.5 ;
1416 sinCPhi = std::sin(cPhi) ;
1417 cosCPhi = std::cos(cPhi) ;
1418
1419 // add angle calculation with correction
1420 // of the difference in domain of atan2 and Sphi
1421 //
1422 vphi = std::atan2(v.y(),v.x()) ;
1423
1424 if ( vphi < fSPhi - kAngTolerance*0.5 ) { vphi += twopi; }
1425 else if ( vphi > fSPhi + fDPhi + kAngTolerance*0.5 ) { vphi -= twopi; }
1426
1427
1428 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1429 {
1430 // pDist -ve when inside
1431
1432 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1433 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1434
1435 // Comp -ve when in direction of outwards normal
1436
1437 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1438 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1439 sidephi = kNull;
1440
1441 if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance)
1442 && (pDistE <= 0.5*kCarTolerance) ) )
1443 || ( (fDPhi > pi) && !((pDistS > 0.5*kCarTolerance)
1444 && (pDistE > 0.5*kCarTolerance) ) ) )
1445 {
1446 // Inside both phi *full* planes
1447
1448 if ( compS < 0 )
1449 {
1450 sphi = pDistS/compS ;
1451
1452 if (sphi >= -0.5*kCarTolerance)
1453 {
1454 xi = p.x() + sphi*v.x() ;
1455 yi = p.y() + sphi*v.y() ;
1456
1457 // Check intersecting with correct half-plane
1458 // (if not -> no intersect)
1459 //
1460 if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance))
1461 { sidephi = kSPhi;
1462 if (((fSPhi-0.5*kAngTolerance)<=vphi)
1463 &&((ePhi+0.5*kAngTolerance)>=vphi))
1464 {
1465 sphi = kInfinity;
1466 }
1467 }
1468 else if ((yi*cosCPhi-xi*sinCPhi)>=0)
1469 {
1470 sphi = kInfinity ;
1471 }
1472 else
1473 {
1474 sidephi = kSPhi ;
1475 if ( pDistS > -kCarTolerance*0.5 )
1476 {
1477 sphi = 0.0 ; // Leave by sphi immediately
1478 }
1479 }
1480 }
1481 else
1482 {
1483 sphi = kInfinity ;
1484 }
1485 }
1486 else
1487 {
1488 sphi = kInfinity ;
1489 }
1490
1491 if ( compE < 0 )
1492 {
1493 sphi2 = pDistE/compE ;
1494
1495 // Only check further if < starting phi intersection
1496 //
1497 if ( (sphi2 > -0.5*kCarTolerance) && (sphi2 < sphi) )
1498 {
1499 xi = p.x() + sphi2*v.x() ;
1500 yi = p.y() + sphi2*v.y() ;
1501
1502 if ((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance))
1503 {
1504 // Leaving via ending phi
1505 //
1506 if(!(((fSPhi-0.5*kAngTolerance)<=vphi)
1507 &&((ePhi+0.5*kAngTolerance)>=vphi)))
1508 {
1509 sidephi = kEPhi ;
1510 if ( pDistE <= -kCarTolerance*0.5 ) { sphi = sphi2 ; }
1511 else { sphi = 0.0 ; }
1512 }
1513 }
1514 else // Check intersecting with correct half-plane
1515
1516 if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1517 {
1518 // Leaving via ending phi
1519 //
1520 sidephi = kEPhi ;
1521 if ( pDistE <= -kCarTolerance*0.5 ) { sphi = sphi2 ; }
1522 else { sphi = 0.0 ; }
1523 }
1524 }
1525 }
1526 }
1527 else
1528 {
1529 sphi = kInfinity ;
1530 }
1531 }
1532 else
1533 {
1534 // On z axis + travel not || to z axis -> if phi of vector direction
1535 // within phi of shape, Step limited by rmax, else Step =0
1536
1537 // vphi = std::atan2(v.y(),v.x()) ;//defined previosly
1538 // G4cout<<"In axis vphi="<<vphi
1539 // <<" Sphi="<<fSPhi<<" Ephi="<<ePhi<<G4endl;
1540 // old if ( (fSPhi < vphi) && (vphi < fSPhi + fDPhi) )
1541 // new : correction for if statement, must be '<='
1542
1543 if ( ((fSPhi-0.5*kAngTolerance) <= vphi)
1544 && (vphi <=( ePhi+0.5*kAngTolerance) ))
1545 {
1546 sphi = kInfinity ;
1547 }
1548 else
1549 {
1550 sidephi = kSPhi ; // arbitrary
1551 sphi = 0.0 ;
1552 }
1553 }
1554 if (sphi < snxt) // Order intersecttions
1555 {
1556 snxt = sphi ;
1557 side = sidephi ;
1558 }
1559 }
1560 if (sr < snxt) // Order intersections
1561 {
1562 snxt = sr ;
1563 side = sider ;
1564 }
1565 }
1566 if (calcNorm)
1567 {
1568 switch(side)
1569 {
1570 case kRMax:
1571 // Note: returned vector not normalised
1572 // (divide by fRMax for unit vector)
1573 //
1574 xi = p.x() + snxt*v.x() ;
1575 yi = p.y() + snxt*v.y() ;
1576 *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1577 *validNorm = true ;
1578 break ;
1579
1580 case kRMin:
1581 *validNorm = false ; // Rmin is inconvex
1582 break ;
1583
1584 case kSPhi:
1585 if ( fDPhi <= pi )
1586 {
1587 *n = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
1588 *validNorm = true ;
1589 }
1590 else
1591 {
1592 *validNorm = false ;
1593 }
1594 break ;
1595
1596 case kEPhi:
1597 if (fDPhi <= pi)
1598 {
1599 *n = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
1600 *validNorm = true ;
1601 }
1602 else
1603 {
1604 *validNorm = false ;
1605 }
1606 break ;
1607
1608 case kPZ:
1609 *n=G4ThreeVector(0,0,1) ;
1610 *validNorm=true ;
1611 break ;
1612
1613 case kMZ:
1614 *n = G4ThreeVector(0,0,-1) ;
1615 *validNorm = true ;
1616 break ;
1617
1618 default:
1619 G4cout.precision(16) ;
1620 G4cout << G4endl ;
1621 DumpInfo();
1622 G4cout << "Position:" << G4endl << G4endl ;
1623 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1624 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1625 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1626 G4cout << "Direction:" << G4endl << G4endl ;
1627 G4cout << "v.x() = " << v.x() << G4endl ;
1628 G4cout << "v.y() = " << v.y() << G4endl ;
1629 G4cout << "v.z() = " << v.z() << G4endl << G4endl ;
1630 G4cout << "Proposed distance :" << G4endl << G4endl ;
1631 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl ;
1632 G4Exception("G4Tubs::DistanceToOut(p,v,..)","Notification",JustWarning,
1633 "Undefined side for valid surface normal to solid.");
1634 break ;
1635 }
1636 }
1637 if ( snxt<kCarTolerance*0.5 ) { snxt=0 ; }
1638
1639 return snxt ;
1640}
1641
1642//////////////////////////////////////////////////////////////////////////
1643//
1644// Calculate distance (<=actual) to closest surface of shape from inside
1645
1646G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const
1647{
1648 G4double safe=0.0, rho, safeR1, safeR2, safeZ ;
1649 G4double safePhi, phiC, cosPhiC, sinPhiC, ePhi ;
1650 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1651
1652#ifdef G4CSGDEBUG
1653 if( Inside(p) == kOutside )
1654 {
1655 G4cout.precision(16) ;
1656 G4cout << G4endl ;
1657 DumpInfo();
1658 G4cout << "Position:" << G4endl << G4endl ;
1659 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1660 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1661 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1662 G4Exception("G4Tubs::DistanceToOut(p)", "Notification", JustWarning,
1663 "Point p is outside !?");
1664 }
1665#endif
1666
1667 if ( fRMin )
1668 {
1669 safeR1 = rho - fRMin ;
1670 safeR2 = fRMax - rho ;
1671
1672 if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1673 else { safe = safeR2 ; }
1674 }
1675 else
1676 {
1677 safe = fRMax - rho ;
1678 }
1679 safeZ = fDz - std::fabs(p.z()) ;
1680
1681 if ( safeZ < safe ) { safe = safeZ ; }
1682
1683 // Check if phi divided, Calc distances closest phi plane
1684 //
1685 if ( fDPhi < twopi )
1686 {
1687 // Above/below central phi of Tubs?
1688
1689 phiC = fSPhi + fDPhi*0.5 ;
1690 cosPhiC = std::cos(phiC) ;
1691 sinPhiC = std::sin(phiC) ;
1692
1693 if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1694 {
1695 safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1696 }
1697 else
1698 {
1699 ePhi = fSPhi + fDPhi ;
1700 safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1701 }
1702 if (safePhi < safe) { safe = safePhi ; }
1703 }
1704 if ( safe < 0 ) { safe = 0 ; }
1705
1706 return safe ;
1707}
1708
1709/////////////////////////////////////////////////////////////////////////
1710//
1711// Create a List containing the transformed vertices
1712// Ordering [0-3] -fDz cross section
1713// [4-7] +fDz cross section such that [0] is below [4],
1714// [1] below [5] etc.
1715// Note:
1716// Caller has deletion resposibility
1717// Potential improvement: For last slice, use actual ending angle
1718// to avoid rounding error problems.
1719
1720G4ThreeVectorList*
1721G4Tubs::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
1722{
1723 G4ThreeVectorList* vertices ;
1724 G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
1725 G4double meshAngle, meshRMax, crossAngle,
1726 cosCrossAngle, sinCrossAngle, sAngle;
1727 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1728 G4int crossSection, noCrossSections;
1729
1730 // Compute no of cross-sections necessary to mesh tube
1731 //
1732 noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
1733
1734 if ( noCrossSections < kMinMeshSections )
1735 {
1736 noCrossSections = kMinMeshSections ;
1737 }
1738 else if (noCrossSections>kMaxMeshSections)
1739 {
1740 noCrossSections = kMaxMeshSections ;
1741 }
1742 // noCrossSections = 4 ;
1743
1744 meshAngle = fDPhi/(noCrossSections - 1) ;
1745 // meshAngle = fDPhi/(noCrossSections) ;
1746
1747 meshRMax = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
1748 meshRMin = fRMin - 100*kCarTolerance ;
1749
1750 // If complete in phi, set start angle such that mesh will be at fRMax
1751 // on the x axis. Will give better extent calculations when not rotated.
1752
1753 if (fDPhi == pi*2.0 && fSPhi == 0 ) { sAngle = -meshAngle*0.5 ; }
1754 else { sAngle = fSPhi ; }
1755
1756 vertices = new G4ThreeVectorList();
1757 vertices->reserve(noCrossSections*4);
1758
1759 if ( vertices )
1760 {
1761 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1762 {
1763 // Compute coordinates of cross section at section crossSection
1764
1765 crossAngle = sAngle + crossSection*meshAngle ;
1766 cosCrossAngle = std::cos(crossAngle) ;
1767 sinCrossAngle = std::sin(crossAngle) ;
1768
1769 rMaxX = meshRMax*cosCrossAngle ;
1770 rMaxY = meshRMax*sinCrossAngle ;
1771
1772 if(meshRMin <= 0.0)
1773 {
1774 rMinX = 0.0 ;
1775 rMinY = 0.0 ;
1776 }
1777 else
1778 {
1779 rMinX = meshRMin*cosCrossAngle ;
1780 rMinY = meshRMin*sinCrossAngle ;
1781 }
1782 vertex0 = G4ThreeVector(rMinX,rMinY,-fDz) ;
1783 vertex1 = G4ThreeVector(rMaxX,rMaxY,-fDz) ;
1784 vertex2 = G4ThreeVector(rMaxX,rMaxY,+fDz) ;
1785 vertex3 = G4ThreeVector(rMinX,rMinY,+fDz) ;
1786
1787 vertices->push_back(pTransform.TransformPoint(vertex0)) ;
1788 vertices->push_back(pTransform.TransformPoint(vertex1)) ;
1789 vertices->push_back(pTransform.TransformPoint(vertex2)) ;
1790 vertices->push_back(pTransform.TransformPoint(vertex3)) ;
1791 }
1792 }
1793 else
1794 {
1795 DumpInfo();
1796 G4Exception("G4Tubs::CreateRotatedVertices()",
1797 "FatalError", FatalException,
1798 "Error in allocation of vertices. Out of memory !");
1799 }
1800 return vertices ;
1801}
1802
1803//////////////////////////////////////////////////////////////////////////
1804//
1805// Stream object contents to an output stream
1806
1807G4GeometryType G4Tubs::GetEntityType() const
1808{
1809 return G4String("G4Tubs");
1810}
1811
1812//////////////////////////////////////////////////////////////////////////
1813//
1814// Stream object contents to an output stream
1815
1816std::ostream& G4Tubs::StreamInfo( std::ostream& os ) const
1817{
1818 os << "-----------------------------------------------------------\n"
1819 << " *** Dump for solid - " << GetName() << " ***\n"
1820 << " ===================================================\n"
1821 << " Solid type: G4Tubs\n"
1822 << " Parameters: \n"
1823 << " inner radius : " << fRMin/mm << " mm \n"
1824 << " outer radius : " << fRMax/mm << " mm \n"
1825 << " half length Z: " << fDz/mm << " mm \n"
1826 << " starting phi : " << fSPhi/degree << " degrees \n"
1827 << " delta phi : " << fDPhi/degree << " degrees \n"
1828 << "-----------------------------------------------------------\n";
1829
1830 return os;
1831}
1832
1833/////////////////////////////////////////////////////////////////////////
1834//
1835// GetPointOnSurface
1836
1837G4ThreeVector G4Tubs::GetPointOnSurface() const
1838{
1839 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1840 aOne, aTwo, aThr, aFou;
1841 G4double rRand;
1842
1843 aOne = 2.*fDz*fDPhi*fRMax;
1844 aTwo = 2.*fDz*fDPhi*fRMin;
1845 aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1846 aFou = 2.*fDz*(fRMax-fRMin);
1847
1848 phi = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1849 cosphi = std::cos(phi);
1850 sinphi = std::sin(phi);
1851
1852 rRand = RandFlat::shoot(fRMin,fRMax);
1853
1854 if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1855
1856 chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1857
1858 if( (chose >=0) && (chose < aOne) )
1859 {
1860 xRand = fRMax*cosphi;
1861 yRand = fRMax*sinphi;
1862 zRand = RandFlat::shoot(-1.*fDz,fDz);
1863 return G4ThreeVector (xRand, yRand, zRand);
1864 }
1865 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1866 {
1867 xRand = fRMin*cosphi;
1868 yRand = fRMin*sinphi;
1869 zRand = RandFlat::shoot(-1.*fDz,fDz);
1870 return G4ThreeVector (xRand, yRand, zRand);
1871 }
1872 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1873 {
1874 xRand = rRand*cosphi;
1875 yRand = rRand*sinphi;
1876 zRand = fDz;
1877 return G4ThreeVector (xRand, yRand, zRand);
1878 }
1879 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1880 {
1881 xRand = rRand*cosphi;
1882 yRand = rRand*sinphi;
1883 zRand = -1.*fDz;
1884 return G4ThreeVector (xRand, yRand, zRand);
1885 }
1886 else if( (chose >= aOne + aTwo + 2.*aThr)
1887 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1888 {
1889 xRand = rRand*std::cos(fSPhi);
1890 yRand = rRand*std::sin(fSPhi);
1891 zRand = RandFlat::shoot(-1.*fDz,fDz);
1892 return G4ThreeVector (xRand, yRand, zRand);
1893 }
1894 else
1895 {
1896 xRand = rRand*std::cos(fSPhi+fDPhi);
1897 yRand = rRand*std::sin(fSPhi+fDPhi);
1898 zRand = RandFlat::shoot(-1.*fDz,fDz);
1899 return G4ThreeVector (xRand, yRand, zRand);
1900 }
1901}
1902
1903///////////////////////////////////////////////////////////////////////////
1904//
1905// Methods for visualisation
1906
1907void G4Tubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1908{
1909 scene.AddSolid (*this) ;
1910}
1911
1912G4Polyhedron* G4Tubs::CreatePolyhedron () const
1913{
1914 return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi) ;
1915}
1916
1917G4NURBS* G4Tubs::CreateNURBS () const
1918{
1919 G4NURBS* pNURBS ;
1920 if (fRMin != 0)
1921 {
1922 if (fDPhi >= twopi)
1923 {
1924 pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ;
1925 }
1926 else
1927 {
1928 pNURBS = new G4NURBStubesector (fRMin,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
1929 }
1930 }
1931 else
1932 {
1933 if (fDPhi >= twopi)
1934 {
1935 pNURBS = new G4NURBScylinder (fRMax,fDz) ;
1936 }
1937 else
1938 {
1939 const G4double epsilon = 1.e-4 ; // Cylinder sector not yet available!
1940 pNURBS = new G4NURBStubesector (epsilon,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
1941 }
1942 }
1943 return pNURBS ;
1944}
Note: See TracBrowser for help on using the repository browser.