| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4Tubs implementation
27 //
28 // 1994-95 P.Kent: first implementation
29 // 08.08.00 V.Grichine: more stable roots of 2
30 // 07.12.00 V.Grichine: phi-section algorithm
31 // 03.05.05 V.Grichine: SurfaceNormal(p) accor
32 // 24.08.16 E.Tcherniaev: reimplemented Calcul
33 // -------------------------------------------
34
35 #include "G4Tubs.hh"
36
37 #if !defined(G4GEOM_USE_UTUBS)
38
39 #include "G4GeomTools.hh"
40 #include "G4VoxelLimits.hh"
41 #include "G4AffineTransform.hh"
42 #include "G4GeometryTolerance.hh"
43 #include "G4BoundingEnvelope.hh"
44
45 #include "G4VPVParameterisation.hh"
46 #include "G4QuickRand.hh"
47
48 #include "G4VGraphicsScene.hh"
49 #include "G4Polyhedron.hh"
50
51 using namespace CLHEP;
52
53 //////////////////////////////////////////////
54 //
55 // Constructor - check parameters, convert ang
56 // - note if pdphi>2PI then reset
57
58 G4Tubs::G4Tubs( const G4String& pName,
59 G4double pRMin, G4double
60 G4double pDz,
61 G4double pSPhi, G4double
62 : G4CSGSolid(pName), fRMin(pRMin), fRMax(pR
63 fSPhi(0), fDPhi(0),
64 fInvRmax( pRMax > 0.0 ? 1.0/pRMax : 0.0 )
65 fInvRmin( pRMin > 0.0 ? 1.0/pRMin : 0.0 )
66 {
67 kRadTolerance = G4GeometryTolerance::GetInst
68 kAngTolerance = G4GeometryTolerance::GetInst
69
70 halfCarTolerance=kCarTolerance*0.5;
71 halfRadTolerance=kRadTolerance*0.5;
72 halfAngTolerance=kAngTolerance*0.5;
73
74 if (pDz<=0) // Check z-len
75 {
76 std::ostringstream message;
77 message << "Negative Z half-length (" << p
78 G4Exception("G4Tubs::G4Tubs()", "GeomSolid
79 }
80 if ( (pRMin >= pRMax) || (pRMin < 0) ) // Ch
81 {
82 std::ostringstream message;
83 message << "Invalid values for radii in so
84 << G4endl
85 << " pRMin = " << pRMin <<
86 G4Exception("G4Tubs::G4Tubs()", "GeomSolid
87 }
88
89 // Check angles
90 //
91 CheckPhiAngles(pSPhi, pDPhi);
92 }
93
94 //////////////////////////////////////////////
95 //
96 // Fake default constructor - sets only member
97 // for usage restri
98 //
99 G4Tubs::G4Tubs( __void__& a )
100 : G4CSGSolid(a)
101 {
102 }
103
104 //////////////////////////////////////////////
105 //
106 // Destructor
107
108 G4Tubs::~G4Tubs() = default;
109
110 //////////////////////////////////////////////
111 //
112 // Copy constructor
113
114 G4Tubs::G4Tubs(const G4Tubs&) = default;
115
116 //////////////////////////////////////////////
117 //
118 // Assignment operator
119
120 G4Tubs& G4Tubs::operator = (const G4Tubs& rhs)
121 {
122 // Check assignment to self
123 //
124 if (this == &rhs) { return *this; }
125
126 // Copy base class data
127 //
128 G4CSGSolid::operator=(rhs);
129
130 // Copy data
131 //
132 kRadTolerance = rhs.kRadTolerance; kAngTole
133 fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz =
134 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
135 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh
136 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r
137 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh
138 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh
139 fPhiFullTube = rhs.fPhiFullTube;
140 fInvRmax = rhs.fInvRmax;
141 fInvRmin = rhs.fInvRmin;
142 halfCarTolerance = rhs.halfCarTolerance;
143 halfRadTolerance = rhs.halfRadTolerance;
144 halfAngTolerance = rhs.halfAngTolerance;
145
146 return *this;
147 }
148
149 //////////////////////////////////////////////
150 //
151 // Dispatch to parameterisation for replicatio
152 // computation & modification.
153
154 void G4Tubs::ComputeDimensions( G4VPVPar
155 const G4int n,
156 const G4VPhysi
157 {
158 p->ComputeDimensions(*this,n,pRep) ;
159 }
160
161 //////////////////////////////////////////////
162 //
163 // Get bounding box
164
165 void G4Tubs::BoundingLimits(G4ThreeVector& pMi
166 {
167 G4double rmin = GetInnerRadius();
168 G4double rmax = GetOuterRadius();
169 G4double dz = GetZHalfLength();
170
171 // Find bounding box
172 //
173 if (GetDeltaPhiAngle() < twopi)
174 {
175 G4TwoVector vmin,vmax;
176 G4GeomTools::DiskExtent(rmin,rmax,
177 GetSinStartPhi(),G
178 GetSinEndPhi(),Get
179 vmin,vmax);
180 pMin.set(vmin.x(),vmin.y(),-dz);
181 pMax.set(vmax.x(),vmax.y(), dz);
182 }
183 else
184 {
185 pMin.set(-rmax,-rmax,-dz);
186 pMax.set( rmax, rmax, dz);
187 }
188
189 // Check correctness of the bounding box
190 //
191 if (pMin.x() >= pMax.x() || pMin.y() >= pMax
192 {
193 std::ostringstream message;
194 message << "Bad bounding box (min >= max)
195 << GetName() << " !"
196 << "\npMin = " << pMin
197 << "\npMax = " << pMax;
198 G4Exception("G4Tubs::BoundingLimits()", "G
199 JustWarning, message);
200 DumpInfo();
201 }
202 }
203
204 //////////////////////////////////////////////
205 //
206 // Calculate extent under transform and specif
207
208 G4bool G4Tubs::CalculateExtent( const EAxis
209 const G4VoxelL
210 const G4Affine
211 G4double
212 G4double
213 {
214 G4ThreeVector bmin, bmax;
215 G4bool exist;
216
217 // Get bounding box
218 BoundingLimits(bmin,bmax);
219
220 // Check bounding box
221 G4BoundingEnvelope bbox(bmin,bmax);
222 #ifdef G4BBOX_EXTENT
223 return bbox.CalculateExtent(pAxis,pVoxelLimi
224 #endif
225 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox
226 {
227 return exist = pMin < pMax;
228 }
229
230 // Get parameters of the solid
231 G4double rmin = GetInnerRadius();
232 G4double rmax = GetOuterRadius();
233 G4double dz = GetZHalfLength();
234 G4double dphi = GetDeltaPhiAngle();
235
236 // Find bounding envelope and calculate exte
237 //
238 const G4int NSTEPS = 24; // numbe
239 G4double astep = twopi/NSTEPS; // max a
240 G4int ksteps = (dphi <= astep) ? 1 : (G4i
241 G4double ang = dphi/ksteps;
242
243 G4double sinHalf = std::sin(0.5*ang);
244 G4double cosHalf = std::cos(0.5*ang);
245 G4double sinStep = 2.*sinHalf*cosHalf;
246 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
247 G4double rext = rmax/cosHalf;
248
249 // bounding envelope for full cylinder consi
250 // in other cases it is a sequence of quadri
251 if (rmin == 0 && dphi == twopi)
252 {
253 G4double sinCur = sinHalf;
254 G4double cosCur = cosHalf;
255
256 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE
257 for (G4int k=0; k<NSTEPS; ++k)
258 {
259 baseA[k].set(rext*cosCur,rext*sinCur,-dz
260 baseB[k].set(rext*cosCur,rext*sinCur, dz
261
262 G4double sinTmp = sinCur;
263 sinCur = sinCur*cosStep + cosCur*sinStep
264 cosCur = cosCur*cosStep - sinTmp*sinStep
265 }
266 std::vector<const G4ThreeVectorList *> pol
267 polygons[0] = &baseA;
268 polygons[1] = &baseB;
269 G4BoundingEnvelope benv(bmin,bmax,polygons
270 exist = benv.CalculateExtent(pAxis,pVoxelL
271 }
272 else
273 {
274 G4double sinStart = GetSinStartPhi();
275 G4double cosStart = GetCosStartPhi();
276 G4double sinEnd = GetSinEndPhi();
277 G4double cosEnd = GetCosEndPhi();
278 G4double sinCur = sinStart*cosHalf + cos
279 G4double cosCur = cosStart*cosHalf - sin
280
281 // set quadrilaterals
282 G4ThreeVectorList pols[NSTEPS+2];
283 for (G4int k=0; k<ksteps+2; ++k) pols[k].r
284 pols[0][0].set(rmin*cosStart,rmin*sinStart
285 pols[0][1].set(rmin*cosStart,rmin*sinStart
286 pols[0][2].set(rmax*cosStart,rmax*sinStart
287 pols[0][3].set(rmax*cosStart,rmax*sinStart
288 for (G4int k=1; k<ksteps+1; ++k)
289 {
290 pols[k][0].set(rmin*cosCur,rmin*sinCur,
291 pols[k][1].set(rmin*cosCur,rmin*sinCur,-
292 pols[k][2].set(rext*cosCur,rext*sinCur,-
293 pols[k][3].set(rext*cosCur,rext*sinCur,
294
295 G4double sinTmp = sinCur;
296 sinCur = sinCur*cosStep + cosCur*sinStep
297 cosCur = cosCur*cosStep - sinTmp*sinStep
298 }
299 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin
300 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin
301 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin
302 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin
303
304 // set envelope and calculate extent
305 std::vector<const G4ThreeVectorList *> pol
306 polygons.resize(ksteps+2);
307 for (G4int k=0; k<ksteps+2; ++k) polygons[
308 G4BoundingEnvelope benv(bmin,bmax,polygons
309 exist = benv.CalculateExtent(pAxis,pVoxelL
310 }
311 return exist;
312 }
313
314 //////////////////////////////////////////////
315 //
316 // Return whether point inside/outside/on surf
317
318 EInside G4Tubs::Inside( const G4ThreeVector& p
319 {
320 G4double r2,pPhi,tolRMin,tolRMax;
321 EInside in = kOutside ;
322
323 if (std::fabs(p.z()) <= fDz - halfCarToleran
324 {
325 r2 = p.x()*p.x() + p.y()*p.y() ;
326
327 if (fRMin != 0.0) { tolRMin = fRMin + half
328 else { tolRMin = 0 ; }
329
330 tolRMax = fRMax - halfRadTolerance ;
331
332 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolR
333 {
334 if ( fPhiFullTube )
335 {
336 in = kInside ;
337 }
338 else
339 {
340 // Try inner tolerant phi boundaries (
341 // if not inside, try outer tolerant p
342
343 if ( (tolRMin==0) && (std::fabs(p.x())
344 && (std::fabs(p.y())
345 {
346 in=kSurface;
347 }
348 else
349 {
350 pPhi = std::atan2(p.y(),p.x()) ;
351 if ( pPhi < -halfAngTolerance ) { p
352
353 if ( fSPhi >= 0 )
354 {
355 if ( (std::fabs(pPhi) < halfAngTol
356 && (std::fabs(fSPhi + fDPhi - tw
357 {
358 pPhi += twopi ; // 0 <= pPhi < 2
359 }
360 if ( (pPhi >= fSPhi + halfAngToler
361 && (pPhi <= fSPhi + fDPhi - half
362 {
363 in = kInside ;
364 }
365 else if ( (pPhi >= fSPhi - halfAng
366 && (pPhi <= fSPhi + fDPhi +
367 {
368 in = kSurface ;
369 }
370 }
371 else // fSPhi < 0
372 {
373 if ( (pPhi <= fSPhi + twopi - half
374 && (pPhi >= fSPhi + fDPhi + hal
375 else if ( (pPhi <= fSPhi + twopi +
376 && (pPhi >= fSPhi + fDPhi
377 {
378 in = kSurface ;
379 }
380 else
381 {
382 in = kInside ;
383 }
384 }
385 }
386 }
387 }
388 else // Try generous boundaries
389 {
390 tolRMin = fRMin - halfRadTolerance ;
391 tolRMax = fRMax + halfRadTolerance ;
392
393 if ( tolRMin < 0 ) { tolRMin = 0; }
394
395 if ( (r2 >= tolRMin*tolRMin) && (r2 <= t
396 {
397 if (fPhiFullTube || (r2 <=halfRadToler
398 { // Continuous
399 in = kSurface ;
400 }
401 else // Try outer tolerant phi boundar
402 {
403 pPhi = std::atan2(p.y(),p.x()) ;
404
405 if ( pPhi < -halfAngTolerance) { pP
406 if ( fSPhi >= 0 )
407 {
408 if ( (std::fabs(pPhi) < halfAngTol
409 && (std::fabs(fSPhi + fDPhi - tw
410 {
411 pPhi += twopi ; // 0 <= pPhi < 2
412 }
413 if ( (pPhi >= fSPhi - halfAngToler
414 && (pPhi <= fSPhi + fDPhi + half
415 {
416 in = kSurface ;
417 }
418 }
419 else // fSPhi < 0
420 {
421 if ( (pPhi <= fSPhi + twopi - half
422 && (pPhi >= fSPhi + fDPhi + half
423 else
424 {
425 in = kSurface ;
426 }
427 }
428 }
429 }
430 }
431 }
432 else if (std::fabs(p.z()) <= fDz + halfCarTo
433 { /
434 r2 = p.x()*p.x() + p.y()*p.y() ;
435 tolRMin = fRMin - halfRadTolerance ;
436 tolRMax = fRMax + halfRadTolerance ;
437
438 if ( tolRMin < 0 ) { tolRMin = 0; }
439
440 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tol
441 {
442 if (fPhiFullTube || (r2 <=halfRadToleran
443 { // Continuous i
444 in = kSurface ;
445 }
446 else // Try outer tolerant phi boundarie
447 {
448 pPhi = std::atan2(p.y(),p.x()) ;
449
450 if ( pPhi < -halfAngTolerance ) { pPh
451 if ( fSPhi >= 0 )
452 {
453 if ( (std::fabs(pPhi) < halfAngToler
454 && (std::fabs(fSPhi + fDPhi - twop
455 {
456 pPhi += twopi ; // 0 <= pPhi < 2pi
457 }
458 if ( (pPhi >= fSPhi - halfAngToleran
459 && (pPhi <= fSPhi + fDPhi + halfAn
460 {
461 in = kSurface;
462 }
463 }
464 else // fSPhi < 0
465 {
466 if ( (pPhi <= fSPhi + twopi - halfAn
467 && (pPhi >= fSPhi + fDPhi + halfA
468 else
469 {
470 in = kSurface ;
471 }
472 }
473 }
474 }
475 }
476 return in;
477 }
478
479 //////////////////////////////////////////////
480 //
481 // Return unit normal of surface closest to p
482 // - note if point on z axis, ignore phi divid
483 // - unsafe if point close to z axis a rmin=0
484
485 G4ThreeVector G4Tubs::SurfaceNormal( const G4T
486 {
487 G4int noSurfaces = 0;
488 G4double rho, pPhi;
489 G4double distZ, distRMin, distRMax;
490 G4double distSPhi = kInfinity, distEPhi = kI
491
492 G4ThreeVector norm, sumnorm(0.,0.,0.);
493 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
494 G4ThreeVector nR, nPs, nPe;
495
496 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
497
498 distRMin = std::fabs(rho - fRMin);
499 distRMax = std::fabs(rho - fRMax);
500 distZ = std::fabs(std::fabs(p.z()) - fDz)
501
502 if (!fPhiFullTube) // Protected against (
503 {
504 if ( rho > halfCarTolerance )
505 {
506 pPhi = std::atan2(p.y(),p.x());
507
508 if (pPhi < fSPhi-halfCarTolerance)
509 else if (pPhi > fSPhi+fDPhi+halfCarToler
510
511 distSPhi = std::fabs( pPhi - fSPhi );
512 distEPhi = std::fabs( pPhi - fSPhi - fDP
513 }
514 else if ( fRMin == 0.0 )
515 {
516 distSPhi = 0.;
517 distEPhi = 0.;
518 }
519 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0
520 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0
521 }
522 if ( rho > halfCarTolerance ) { nR = G4Three
523
524 if( distRMax <= halfCarTolerance )
525 {
526 ++noSurfaces;
527 sumnorm += nR;
528 }
529 if( (fRMin != 0.0) && (distRMin <= halfCarTo
530 {
531 ++noSurfaces;
532 sumnorm -= nR;
533 }
534 if( fDPhi < twopi )
535 {
536 if (distSPhi <= halfAngTolerance)
537 {
538 ++noSurfaces;
539 sumnorm += nPs;
540 }
541 if (distEPhi <= halfAngTolerance)
542 {
543 ++noSurfaces;
544 sumnorm += nPe;
545 }
546 }
547 if (distZ <= halfCarTolerance)
548 {
549 ++noSurfaces;
550 if ( p.z() >= 0.) { sumnorm += nZ; }
551 else { sumnorm -= nZ; }
552 }
553 if ( noSurfaces == 0 )
554 {
555 #ifdef G4CSGDEBUG
556 G4Exception("G4Tubs::SurfaceNormal(p)", "G
557 JustWarning, "Point p is not o
558 G4long oldprc = G4cout.precision(20);
559 G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y
560 << G4endl << G4endl;
561 G4cout.precision(oldprc) ;
562 #endif
563 norm = ApproxSurfaceNormal(p);
564 }
565 else if ( noSurfaces == 1 ) { norm = sumnor
566 else { norm = sumnor
567
568 return norm;
569 }
570
571 //////////////////////////////////////////////
572 //
573 // Algorithm for SurfaceNormal() following the
574 // for points not on the surface
575
576 G4ThreeVector G4Tubs::ApproxSurfaceNormal( con
577 {
578 ENorm side ;
579 G4ThreeVector norm ;
580 G4double rho, phi ;
581 G4double distZ, distRMin, distRMax, distSPhi
582
583 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
584
585 distRMin = std::fabs(rho - fRMin) ;
586 distRMax = std::fabs(rho - fRMax) ;
587 distZ = std::fabs(std::fabs(p.z()) - fDz)
588
589 if (distRMin < distRMax) // First minimum
590 {
591 if ( distZ < distRMin )
592 {
593 distMin = distZ ;
594 side = kNZ ;
595 }
596 else
597 {
598 distMin = distRMin ;
599 side = kNRMin ;
600 }
601 }
602 else
603 {
604 if ( distZ < distRMax )
605 {
606 distMin = distZ ;
607 side = kNZ ;
608 }
609 else
610 {
611 distMin = distRMax ;
612 side = kNRMax ;
613 }
614 }
615 if (!fPhiFullTube && (rho != 0.0) ) // Pro
616 {
617 phi = std::atan2(p.y(),p.x()) ;
618
619 if ( phi < 0 ) { phi += twopi; }
620
621 if ( fSPhi < 0 )
622 {
623 distSPhi = std::fabs(phi - (fSPhi + twop
624 }
625 else
626 {
627 distSPhi = std::fabs(phi - fSPhi)*rho ;
628 }
629 distEPhi = std::fabs(phi - fSPhi - fDPhi)*
630
631 if (distSPhi < distEPhi) // Find new minim
632 {
633 if ( distSPhi < distMin )
634 {
635 side = kNSPhi ;
636 }
637 }
638 else
639 {
640 if ( distEPhi < distMin )
641 {
642 side = kNEPhi ;
643 }
644 }
645 }
646 switch ( side )
647 {
648 case kNRMin : // Inner radius
649 {
650 norm = G4ThreeVector(-p.x()/rho, -p.y()/
651 break ;
652 }
653 case kNRMax : // Outer radius
654 {
655 norm = G4ThreeVector(p.x()/rho, p.y()/rh
656 break ;
657 }
658 case kNZ : // + or - dz
659 {
660 if ( p.z() > 0 ) { norm = G4ThreeVector
661 else { norm = G4ThreeVector
662 break ;
663 }
664 case kNSPhi:
665 {
666 norm = G4ThreeVector(sinSPhi, -cosSPhi,
667 break ;
668 }
669 case kNEPhi:
670 {
671 norm = G4ThreeVector(-sinEPhi, cosEPhi,
672 break;
673 }
674 default: // Should never reach this c
675 {
676 DumpInfo();
677 G4Exception("G4Tubs::ApproxSurfaceNormal
678 "GeomSolids1002", JustWarnin
679 "Undefined side for valid su
680 break ;
681 }
682 }
683 return norm;
684 }
685
686 //////////////////////////////////////////////
687 //
688 //
689 // Calculate distance to shape from outside, a
690 // - return kInfinity if no intersection, or i
691 //
692 // - Compute the intersection with the z plane
693 // - if at valid r, phi, return
694 //
695 // -> If point is outer outer radius, compute
696 // - if at valid phi,z return
697 //
698 // -> Compute intersection with inner radius,
699 // - if valid (in z,phi), save intersct
700 //
701 // -> If phi segmented, compute intersectio
702 // - return smallest of valid phi inter
703 // inner radius intersection
704 //
705 // NOTE:
706 // - 'if valid' implies tolerant checking of i
707
708 G4double G4Tubs::DistanceToIn( const G4ThreeVe
709 const G4ThreeVe
710 {
711 G4double snxt = kInfinity ; // snxt = d
712 G4double tolORMin2, tolIRMax2 ; // 'generou
713 G4double tolORMax2, tolIRMin2, tolODz, tolID
714 const G4double dRmax = 100.*fRMax;
715
716 // Intersection point variables
717 //
718 G4double Dist, sd, xi, yi, zi, rho2, inum, i
719 G4double t1, t2, t3, b, c, d ; // Quadra
720
721 // Calculate tolerant rmin and rmax
722
723 if (fRMin > kRadTolerance)
724 {
725 tolORMin2 = (fRMin - halfRadTolerance)*(fR
726 tolIRMin2 = (fRMin + halfRadTolerance)*(fR
727 }
728 else
729 {
730 tolORMin2 = 0.0 ;
731 tolIRMin2 = 0.0 ;
732 }
733 tolORMax2 = (fRMax + halfRadTolerance)*(fRMa
734 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMa
735
736 // Intersection with Z surfaces
737
738 tolIDz = fDz - halfCarTolerance ;
739 tolODz = fDz + halfCarTolerance ;
740
741 if (std::fabs(p.z()) >= tolIDz)
742 {
743 if ( p.z()*v.z() < 0 ) // at +Z going i
744 {
745 sd = (std::fabs(p.z()) - fDz)/std::fabs(
746
747 if(sd < 0.0) { sd = 0.0; }
748
749 xi = p.x() + sd*v.x() ;
750 yi = p.y() + sd*v.y() ;
751 rho2 = xi*xi + yi*yi ;
752
753 // Check validity of intersection
754
755 if ((tolIRMin2 <= rho2) && (rho2 <= tolI
756 {
757 if (!fPhiFullTube && (rho2 != 0.0))
758 {
759 // Psi = angle made with central (av
760 //
761 inum = xi*cosCPhi + yi*sinCPhi ;
762 iden = std::sqrt(rho2) ;
763 cosPsi = inum/iden ;
764 if (cosPsi >= cosHDPhiIT) { return
765 }
766 else
767 {
768 return sd ;
769 }
770 }
771 }
772 else
773 {
774 if ( snxt<halfCarTolerance ) { snxt=0;
775 return snxt ; // On/outside extent, and
776 // -> cannot intersect
777 }
778 }
779
780 // -> Can not intersect z surfaces
781 //
782 // Intersection with rmax (possible return)
783 //
784 // Intersection point (xi,yi,zi) on line x=p
785 //
786 // Intersects with x^2+y^2=R^2
787 //
788 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.
789 // t1 t2
790
791 t1 = 1.0 - v.z()*v.z() ;
792 t2 = p.x()*v.x() + p.y()*v.y() ;
793 t3 = p.x()*p.x() + p.y()*p.y() ;
794
795 if ( t1 > 0 ) // Check not || to z ax
796 {
797 b = t2/t1 ;
798 c = t3 - fRMax*fRMax ;
799 if ((t3 >= tolORMax2) && (t2<0)) // This
800 {
801 // Try outer cylinder intersection
802 // c=(t3-fRMax*fRMax)/t1;
803
804 c /= t1 ;
805 d = b*b - c ;
806
807 if (d >= 0) // If real root
808 {
809 sd = c/(-b+std::sqrt(d));
810 if (sd >= 0) // If 'forwards'
811 {
812 if ( sd>dRmax ) // Avoid rounding er
813 { // 64 bits systems.
814 G4double fTerm = sd-std::fmod(sd,d
815 sd = fTerm + DistanceToIn(p+fTerm*
816 }
817 // Check z intersection
818 //
819 zi = p.z() + sd*v.z() ;
820 if (std::fabs(zi)<=tolODz)
821 {
822 // Z ok. Check phi intersection if
823 //
824 if (fPhiFullTube)
825 {
826 return sd ;
827 }
828 else
829 {
830 xi = p.x() + sd*v.x() ;
831 yi = p.y() + sd*v.y() ;
832 cosPsi = (xi*cosCPhi + yi*sinCPh
833 if (cosPsi >= cosHDPhiIT) { ret
834 }
835 } // end if std::fabs(zi)
836 } // end if (sd>=0)
837 } // end if (d>=0)
838 } // end if (r>=fRMax)
839 else
840 {
841 // Inside outer radius :
842 // check not inside, and heading through
843
844 if ((t3 > tolIRMin2) && (t2 < 0) && (std
845 {
846 // Inside both radii, delta r -ve, ins
847
848 if (!fPhiFullTube)
849 {
850 inum = p.x()*cosCPhi + p.y()*sinCP
851 iden = std::sqrt(t3) ;
852 cosPsi = inum/iden ;
853 if (cosPsi >= cosHDPhiIT)
854 {
855 // In the old version, the small n
856 // on surface was not taken in acc
857 // New version: check the tangent
858 // if no intersection, return kInf
859 // return sd.
860 //
861 c = t3-fRMax*fRMax;
862 if ( c<=0.0 )
863 {
864 return 0.0;
865 }
866 else
867 {
868 c = c/t1 ;
869 d = b*b-c;
870 if ( d>=0.0 )
871 {
872 snxt = c/(-b+std::sqrt(d)); //
873 //
874 if ( snxt < halfCarTolerance )
875 return snxt ;
876 }
877 else
878 {
879 return kInfinity;
880 }
881 }
882 }
883 }
884 else
885 {
886 // In the old version, the small neg
887 // on surface was not taken in accou
888 // New version: check the tangent fo
889 // if no intersection, return kInfin
890 // return sd.
891 //
892 c = t3 - fRMax*fRMax;
893 if ( c<=0.0 )
894 {
895 return 0.0;
896 }
897 else
898 {
899 c = c/t1 ;
900 d = b*b-c;
901 if ( d>=0.0 )
902 {
903 snxt= c/(-b+std::sqrt(d)); // us
904 // fo
905 if ( snxt < halfCarTolerance ) {
906 return snxt ;
907 }
908 else
909 {
910 return kInfinity;
911 }
912 }
913 } // end if (!fPhiFullTube)
914 } // end if (t3>tolIRMin2)
915 } // end if (Inside Outer Radius)
916 if ( fRMin != 0.0 ) // Try inner cylind
917 {
918 c = (t3 - fRMin*fRMin)/t1 ;
919 d = b*b - c ;
920 if ( d >= 0.0 ) // If real root
921 {
922 // Always want 2nd root - we are outsi
923 // - If on surface of rmin also need f
924
925 sd =( b > 0. )? c/(-b - std::sqrt(d))
926 if (sd >= -halfCarTolerance) // check
927 {
928 // Check z intersection
929 //
930 if(sd < 0.0) { sd = 0.0; }
931 if ( sd>dRmax ) // Avoid rounding er
932 { // 64 bits systems.
933 G4double fTerm = sd-std::fmod(sd,d
934 sd = fTerm + DistanceToIn(p+fTerm*
935 }
936 zi = p.z() + sd*v.z() ;
937 if (std::fabs(zi) <= tolODz)
938 {
939 // Z ok. Check phi
940 //
941 if ( fPhiFullTube )
942 {
943 return sd ;
944 }
945 else
946 {
947 xi = p.x() + sd*v.x() ;
948 yi = p.y() + sd*v.y() ;
949 cosPsi = (xi*cosCPhi + yi*sinCPh
950 if (cosPsi >= cosHDPhiIT)
951 {
952 // Good inner radius isect
953 // - but earlier phi isect sti
954
955 snxt = sd ;
956 }
957 }
958 } // end if std::fabs(zi)
959 } // end if (sd>=0)
960 } // end if (d>=0)
961 } // end if (fRMin)
962 }
963
964 // Phi segment intersection
965 //
966 // o Tolerant of points inside phi planes by
967 //
968 // o NOTE: Large duplication of code between
969 // -> only diffs: sphi -> ephi, Comp
970 // intersection check <=0 -> >=0
971 // -> use some form of loop Construc
972 //
973 if ( !fPhiFullTube )
974 {
975 // First phi surface (Starting phi)
976 //
977 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
978
979 if ( Comp < 0 ) // Component in outwards
980 {
981 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
982
983 if ( Dist < halfCarTolerance )
984 {
985 sd = Dist/Comp ;
986
987 if (sd < snxt)
988 {
989 if ( sd < 0 ) { sd = 0.0; }
990 zi = p.z() + sd*v.z() ;
991 if ( std::fabs(zi) <= tolODz )
992 {
993 xi = p.x() + sd*v.x() ;
994 yi = p.y() + sd*v.y() ;
995 rho2 = xi*xi + yi*yi ;
996
997 if ( ( (rho2 >= tolIRMin2) && (rho
998 || ( (rho2 > tolORMin2) && (rho
999 && ( v.y()*cosSPhi - v.x()*sin
1000 && ( v.x()*cosSPhi + v.y()*si
1001 || ( (rho2 > tolIRMax2) && (rho
1002 && (v.y()*cosSPhi - v.x()*sin
1003 && (v.x()*cosSPhi + v.y()*sin
1004 {
1005 // z and r intersections good
1006 // - check intersecting with co
1007 //
1008 if ((yi*cosCPhi-xi*sinCPhi) <=
1009 }
1010 }
1011 }
1012 }
1013 }
1014
1015 // Second phi surface (Ending phi)
1016
1017 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi
1018
1019 if (Comp < 0 ) // Component in outwards
1020 {
1021 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi)
1022
1023 if ( Dist < halfCarTolerance )
1024 {
1025 sd = Dist/Comp ;
1026
1027 if (sd < snxt)
1028 {
1029 if ( sd < 0 ) { sd = 0; }
1030 zi = p.z() + sd*v.z() ;
1031 if ( std::fabs(zi) <= tolODz )
1032 {
1033 xi = p.x() + sd*v.x() ;
1034 yi = p.y() + sd*v.y() ;
1035 rho2 = xi*xi + yi*yi ;
1036 if ( ( (rho2 >= tolIRMin2) && (rh
1037 || ( (rho2 > tolORMin2) && (
1038 && (v.x()*sinEPhi - v.y()*c
1039 && (v.x()*cosEPhi + v.y()*s
1040 || ( (rho2 > tolIRMax2) && (r
1041 && (v.x()*sinEPhi - v.y()*c
1042 && (v.x()*cosEPhi + v.y()*s
1043 {
1044 // z and r intersections good
1045 // - check intersecting with co
1046 //
1047 if ( (yi*cosCPhi-xi*sinCPhi) >=
1048 } //?? >=
1049 }
1050 }
1051 }
1052 } // Comp < 0
1053 } // !fPhiFullTube
1054 if ( snxt<halfCarTolerance ) { snxt=0; }
1055 return snxt ;
1056 }
1057
1058 /////////////////////////////////////////////
1059 //
1060 // Calculate distance to shape from outside,
1061 // - return kInfinity if no intersection, or
1062 //
1063 // - Compute the intersection with the z plan
1064 // - if at valid r, phi, return
1065 //
1066 // -> If point is outer outer radius, compute
1067 // - if at valid phi,z return
1068 //
1069 // -> Compute intersection with inner radius,
1070 // - if valid (in z,phi), save intersc
1071 //
1072 // -> If phi segmented, compute intersecti
1073 // - return smallest of valid phi inte
1074 // inner radius intersection
1075 //
1076 // NOTE:
1077 // - Precalculations for phi trigonometry are
1078 // - `if valid' implies tolerant checking of
1079 // Calculate distance (<= actual) to closes
1080 // - Calculate distance to z, radial planes
1081 // - Only to phi planes if outside phi extent
1082 // - Return 0 if point inside
1083
1084 G4double G4Tubs::DistanceToIn( const G4ThreeV
1085 {
1086 G4double safe=0.0, rho, safe1, safe2, safe3
1087 G4double safePhi, cosPsi ;
1088
1089 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()
1090 safe1 = fRMin - rho ;
1091 safe2 = rho - fRMax ;
1092 safe3 = std::fabs(p.z()) - fDz ;
1093
1094 if ( safe1 > safe2 ) { safe = safe1; }
1095 else { safe = safe2; }
1096 if ( safe3 > safe ) { safe = safe3; }
1097
1098 if ( (!fPhiFullTube) && ((rho) != 0.0) )
1099 {
1100 // Psi=angle from central phi to point
1101 //
1102 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/
1103
1104 if ( cosPsi < cosHDPhi )
1105 {
1106 // Point lies outside phi range
1107
1108 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <=
1109 {
1110 safePhi = std::fabs(p.x()*sinSPhi - p
1111 }
1112 else
1113 {
1114 safePhi = std::fabs(p.x()*sinEPhi - p
1115 }
1116 if ( safePhi > safe ) { safe = safePhi
1117 }
1118 }
1119 if ( safe < 0 ) { safe = 0; }
1120 return safe ;
1121 }
1122
1123 /////////////////////////////////////////////
1124 //
1125 // Calculate distance to surface of shape fro
1126 // - Only Calc rmax intersection if no valid
1127
1128 G4double G4Tubs::DistanceToOut( const G4Three
1129 const G4Three
1130 const G4bool
1131 G4bool*
1132 G4Three
1133 {
1134 ESide side=kNull , sider=kNull, sidephi=kNu
1135 G4double snxt, srd=kInfinity, sphi=kInfinit
1136 G4double deltaR, t1, t2, t3, b, c, d2, roMi
1137
1138 // Vars for phi intersection:
1139
1140 G4double pDistS, compS, pDistE, compE, sphi
1141
1142 // Z plane intersection
1143
1144 if (v.z() > 0 )
1145 {
1146 pdist = fDz - p.z() ;
1147 if ( pdist > halfCarTolerance )
1148 {
1149 snxt = pdist/v.z() ;
1150 side = kPZ ;
1151 }
1152 else
1153 {
1154 if (calcNorm)
1155 {
1156 *n = G4ThreeVector(0,0,1) ;
1157 *validNorm = true ;
1158 }
1159 return snxt = 0 ;
1160 }
1161 }
1162 else if ( v.z() < 0 )
1163 {
1164 pdist = fDz + p.z() ;
1165
1166 if ( pdist > halfCarTolerance )
1167 {
1168 snxt = -pdist/v.z() ;
1169 side = kMZ ;
1170 }
1171 else
1172 {
1173 if (calcNorm)
1174 {
1175 *n = G4ThreeVector(0,0,-1) ;
1176 *validNorm = true ;
1177 }
1178 return snxt = 0.0 ;
1179 }
1180 }
1181 else
1182 {
1183 snxt = kInfinity ; // Travel perpendic
1184 side = kNull;
1185 }
1186
1187 // Radial Intersections
1188 //
1189 // Find intersection with cylinders at rmax
1190 // Intersection point (xi,yi,zi) on line x=
1191 //
1192 // Intersects with x^2+y^2=R^2
1193 //
1194 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v
1195 //
1196 // t1 t2
1197
1198 t1 = 1.0 - v.z()*v.z() ; // since v
1199 t2 = p.x()*v.x() + p.y()*v.y() ;
1200 t3 = p.x()*p.x() + p.y()*p.y() ;
1201
1202 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fR
1203 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t
1204
1205 if ( t1 > 0 ) // Check not parallel
1206 {
1207 // Calculate srd, r exit distance
1208
1209 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax
1210 {
1211 // Delta r not negative => leaving via
1212
1213 deltaR = t3 - fRMax*fRMax ;
1214
1215 // NOTE: Should use rho-fRMax<-kRadTole
1216 // - avoid sqrt for efficiency
1217
1218 if ( deltaR < -kRadTolerance*fRMax )
1219 {
1220 b = t2/t1 ;
1221 c = deltaR/t1 ;
1222 d2 = b*b-c;
1223 if( d2 >= 0 ) { srd = c/( -b - std::s
1224 else { srd = 0.; }
1225 sider = kRMax ;
1226 }
1227 else
1228 {
1229 // On tolerant boundary & heading out
1230 // outer radial surface -> leaving im
1231
1232 if ( calcNorm )
1233 {
1234 G4double invRho = FastInverseRxy( p
1235 *n = G4ThreeVector(p.x()*in
1236 *validNorm = true ;
1237 }
1238 return snxt = 0 ; // Leaving by rmax
1239 }
1240 }
1241 else if ( t2 < 0. ) // i.e. t2 < 0; Poss
1242 {
1243 roMin2 = t3 - t2*t2/t1 ; // min ro2 of
1244
1245 if ( (fRMin != 0.0) && (roMin2 < fRMin*
1246 {
1247 deltaR = t3 - fRMin*fRMin ;
1248 b = t2/t1 ;
1249 c = deltaR/t1 ;
1250 d2 = b*b - c ;
1251
1252 if ( d2 >= 0 ) // Leaving via rmin
1253 {
1254 // NOTE: SHould use rho-rmin>kRadTo
1255 // - avoid sqrt for efficiency
1256
1257 if (deltaR > kRadTolerance*fRMin)
1258 {
1259 srd = c/(-b+std::sqrt(d2));
1260 sider = kRMin ;
1261 }
1262 else
1263 {
1264 if ( calcNorm ) {
1265 *validNorm = false;
1266 } // Concave side
1267 return snxt = 0.0;
1268 }
1269 }
1270 else // No rmin intersect -> must
1271 {
1272 deltaR = t3 - fRMax*fRMax ;
1273 c = deltaR/t1 ;
1274 d2 = b*b-c;
1275 if( d2 >=0. )
1276 {
1277 srd = -b + std::sqrt(d2) ;
1278 sider = kRMax ;
1279 }
1280 else // Case: On the border+t2<kRad
1281 // (v is perpendicular t
1282 {
1283 if (calcNorm)
1284 {
1285 G4double invRho = FastInverseRx
1286 *n = G4ThreeVector(p.x()*invRho
1287 *validNorm = true ;
1288 }
1289 return snxt = 0.0;
1290 }
1291 }
1292 }
1293 else if ( roi2 > fRMax*(fRMax + kRadTol
1294 // No rmin intersect -> must be rm
1295 {
1296 deltaR = t3 - fRMax*fRMax ;
1297 b = t2/t1 ;
1298 c = deltaR/t1;
1299 d2 = b*b-c;
1300 if( d2 >= 0 )
1301 {
1302 srd = -b + std::sqrt(d2) ;
1303 sider = kRMax ;
1304 }
1305 else // Case: On the border+t2<kRadTo
1306 // (v is perpendicular to
1307 {
1308 if (calcNorm)
1309 {
1310 G4double invRho = FastInverseRxy(
1311 *n = G4ThreeVector(p.x()*invRho,p
1312 *validNorm = true ;
1313 }
1314 return snxt = 0.0;
1315 }
1316 }
1317 }
1318
1319 // Phi Intersection
1320
1321 if ( !fPhiFullTube )
1322 {
1323 // add angle calculation with correctio
1324 // of the difference in domain of atan2
1325 //
1326 vphi = std::atan2(v.y(),v.x()) ;
1327
1328 if ( vphi < fSPhi - halfAngTolerance )
1329 else if ( vphi > fSPhi + fDPhi + halfAn
1330
1331
1332 if ( (p.x() != 0.0) || (p.y() != 0.0) )
1333 {
1334 // pDist -ve when inside
1335
1336 pDistS = p.x()*sinSPhi - p.y()*cosSPh
1337 pDistE = -p.x()*sinEPhi + p.y()*cosEP
1338
1339 // Comp -ve when in direction of outw
1340
1341 compS = -sinSPhi*v.x() + cosSPhi*v.y(
1342 compE = sinEPhi*v.x() - cosEPhi*v.y(
1343
1344 sidephi = kNull;
1345
1346 if( ( (fDPhi <= pi) && ( (pDistS <= h
1347 && (pDistE <= h
1348 || ( (fDPhi > pi) && ((pDistS <= h
1349 || (pDistE <=
1350 {
1351 // Inside both phi *full* planes
1352
1353 if ( compS < 0 )
1354 {
1355 sphi = pDistS/compS ;
1356
1357 if (sphi >= -halfCarTolerance)
1358 {
1359 xi = p.x() + sphi*v.x() ;
1360 yi = p.y() + sphi*v.y() ;
1361
1362 // Check intersecting with corr
1363 // (if not -> no intersect)
1364 //
1365 if((std::fabs(xi)<=kCarToleranc
1366 {
1367 sidephi = kSPhi;
1368 if (((fSPhi-halfAngTolerance)
1369 &&((fSPhi+fDPhi+halfAngTol
1370 {
1371 sphi = kInfinity;
1372 }
1373 }
1374 else if ( yi*cosCPhi-xi*sinCPhi
1375 {
1376 sphi = kInfinity ;
1377 }
1378 else
1379 {
1380 sidephi = kSPhi ;
1381 if ( pDistS > -halfCarToleran
1382 {
1383 sphi = 0.0 ; // Leave by sp
1384 }
1385 }
1386 }
1387 else
1388 {
1389 sphi = kInfinity ;
1390 }
1391 }
1392 else
1393 {
1394 sphi = kInfinity ;
1395 }
1396
1397 if ( compE < 0 )
1398 {
1399 sphi2 = pDistE/compE ;
1400
1401 // Only check further if < starti
1402 //
1403 if ( (sphi2 > -halfCarTolerance)
1404 {
1405 xi = p.x() + sphi2*v.x() ;
1406 yi = p.y() + sphi2*v.y() ;
1407
1408 if((std::fabs(xi)<=kCarToleranc
1409 {
1410 // Leaving via ending phi
1411 //
1412 if( (fSPhi-halfAngTolerance >
1413 ||(fSPhi+fDPhi+halfAngTo
1414 {
1415 sidephi = kEPhi ;
1416 if ( pDistE <= -halfCarTole
1417 else
1418 }
1419 }
1420 else // Check intersecting w
1421
1422 if ( (yi*cosCPhi-xi*sinCPhi) >=
1423 {
1424 // Leaving via ending phi
1425 //
1426 sidephi = kEPhi ;
1427 if ( pDistE <= -halfCarTolera
1428 else
1429 }
1430 }
1431 }
1432 }
1433 else
1434 {
1435 sphi = kInfinity ;
1436 }
1437 }
1438 else
1439 {
1440 // On z axis + travel not || to z axi
1441 // within phi of shape, Step limited
1442
1443 if ( (fSPhi - halfAngTolerance <= vph
1444 && (vphi <= fSPhi + fDPhi + halfAn
1445 {
1446 sphi = kInfinity ;
1447 }
1448 else
1449 {
1450 sidephi = kSPhi ; // arbitrary
1451 sphi = 0.0 ;
1452 }
1453 }
1454 if (sphi < snxt) // Order intersecttio
1455 {
1456 snxt = sphi ;
1457 side = sidephi ;
1458 }
1459 }
1460 if (srd < snxt) // Order intersections
1461 {
1462 snxt = srd ;
1463 side = sider ;
1464 }
1465 }
1466 if (calcNorm)
1467 {
1468 switch(side)
1469 {
1470 case kRMax:
1471 // Note: returned vector not normalis
1472 // (divide by fRMax for unit vector)
1473 //
1474 xi = p.x() + snxt*v.x() ;
1475 yi = p.y() + snxt*v.y() ;
1476 *n = G4ThreeVector(xi/fRMax,yi/fRMax,
1477 *validNorm = true ;
1478 break ;
1479
1480 case kRMin:
1481 *validNorm = false ; // Rmin is inco
1482 break ;
1483
1484 case kSPhi:
1485 if ( fDPhi <= pi )
1486 {
1487 *n = G4ThreeVector(sinSPhi,
1488 *validNorm = true ;
1489 }
1490 else
1491 {
1492 *validNorm = false ;
1493 }
1494 break ;
1495
1496 case kEPhi:
1497 if (fDPhi <= pi)
1498 {
1499 *n = G4ThreeVector(-sinEPhi,cosEPhi
1500 *validNorm = true ;
1501 }
1502 else
1503 {
1504 *validNorm = false ;
1505 }
1506 break ;
1507
1508 case kPZ:
1509 *n = G4ThreeVector(0,0,1) ;
1510 *validNorm = true ;
1511 break ;
1512
1513 case kMZ:
1514 *n = G4ThreeVector(0,0,-1) ;
1515 *validNorm = true ;
1516 break ;
1517
1518 default:
1519 G4cout << G4endl ;
1520 DumpInfo();
1521 std::ostringstream message;
1522 G4long oldprc = message.precision(16)
1523 message << "Undefined side for valid
1524 << G4endl
1525 << "Position:" << G4endl <<
1526 << "p.x() = " << p.x()/mm <
1527 << "p.y() = " << p.y()/mm <
1528 << "p.z() = " << p.z()/mm <
1529 << "Direction:" << G4endl <<
1530 << "v.x() = " << v.x() << G
1531 << "v.y() = " << v.y() << G
1532 << "v.z() = " << v.z() << G
1533 << "Proposed distance :" << G
1534 << "snxt = " << snxt/mm <<
1535 message.precision(oldprc) ;
1536 G4Exception("G4Tubs::DistanceToOut(p,
1537 JustWarning, message);
1538 break ;
1539 }
1540 }
1541 if ( snxt<halfCarTolerance ) { snxt=0 ; }
1542
1543 return snxt ;
1544 }
1545
1546 /////////////////////////////////////////////
1547 //
1548 // Calculate distance (<=actual) to closest s
1549
1550 G4double G4Tubs::DistanceToOut( const G4Three
1551 {
1552 G4double safe=0.0, rho, safeR1, safeR2, saf
1553 rho = std::sqrt(p.x()*p.x() + p.y()*p.y())
1554
1555 #ifdef G4CSGDEBUG
1556 if( Inside(p) == kOutside )
1557 {
1558 G4long oldprc = G4cout.precision(16) ;
1559 G4cout << G4endl ;
1560 DumpInfo();
1561 G4cout << "Position:" << G4endl << G4end
1562 G4cout << "p.x() = " << p.x()/mm << " m
1563 G4cout << "p.y() = " << p.y()/mm << " m
1564 G4cout << "p.z() = " << p.z()/mm << " m
1565 G4cout.precision(oldprc) ;
1566 G4Exception("G4Tubs::DistanceToOut(p)", "
1567 JustWarning, "Point p is outs
1568 }
1569 #endif
1570
1571 if ( fRMin != 0.0 )
1572 {
1573 safeR1 = rho - fRMin ;
1574 safeR2 = fRMax - rho ;
1575
1576 if ( safeR1 < safeR2 ) { safe = safeR1 ;
1577 else { safe = safeR2 ;
1578 }
1579 else
1580 {
1581 safe = fRMax - rho ;
1582 }
1583 safeZ = fDz - std::fabs(p.z()) ;
1584
1585 if ( safeZ < safe ) { safe = safeZ ; }
1586
1587 // Check if phi divided, Calc distances clo
1588 //
1589 if ( !fPhiFullTube )
1590 {
1591 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1592 {
1593 safePhi = -(p.x()*sinSPhi - p.y()*cosSP
1594 }
1595 else
1596 {
1597 safePhi = (p.x()*sinEPhi - p.y()*cosEPh
1598 }
1599 if (safePhi < safe) { safe = safePhi ; }
1600 }
1601 if ( safe < 0 ) { safe = 0 ; }
1602
1603 return safe ;
1604 }
1605
1606 /////////////////////////////////////////////
1607 //
1608 // Stream object contents to an output stream
1609
1610 G4GeometryType G4Tubs::GetEntityType() const
1611 {
1612 return {"G4Tubs"};
1613 }
1614
1615 /////////////////////////////////////////////
1616 //
1617 // Make a clone of the object
1618 //
1619 G4VSolid* G4Tubs::Clone() const
1620 {
1621 return new G4Tubs(*this);
1622 }
1623
1624 /////////////////////////////////////////////
1625 //
1626 // Stream object contents to an output stream
1627
1628 std::ostream& G4Tubs::StreamInfo( std::ostrea
1629 {
1630 G4long oldprc = os.precision(16);
1631 os << "------------------------------------
1632 << " *** Dump for solid - " << GetNam
1633 << " ================================
1634 << " Solid type: G4Tubs\n"
1635 << " Parameters: \n"
1636 << " inner radius : " << fRMin/mm <<
1637 << " outer radius : " << fRMax/mm <<
1638 << " half length Z: " << fDz/mm << "
1639 << " starting phi : " << fSPhi/degree
1640 << " delta phi : " << fDPhi/degree
1641 << "------------------------------------
1642 os.precision(oldprc);
1643
1644 return os;
1645 }
1646
1647 /////////////////////////////////////////////
1648 //
1649 // GetPointOnSurface
1650
1651 G4ThreeVector G4Tubs::GetPointOnSurface() con
1652 {
1653 G4double Rmax = fRMax;
1654 G4double Rmin = fRMin;
1655 G4double hz = 2.*fDz; // height
1656 G4double lext = fDPhi*Rmax; // length of ex
1657 G4double lint = fDPhi*Rmin; // length of in
1658
1659 // Set array of surface areas
1660 //
1661 G4double RRmax = Rmax * Rmax;
1662 G4double RRmin = Rmin * Rmin;
1663 G4double sbase = 0.5*fDPhi*(RRmax - RRmin);
1664 G4double scut = (fDPhi == twopi) ? 0. : hz*
1665 G4double ssurf[6] = { scut, scut, sbase, sb
1666 ssurf[1] += ssurf[0];
1667 ssurf[2] += ssurf[1];
1668 ssurf[3] += ssurf[2];
1669 ssurf[4] += ssurf[3];
1670 ssurf[5] += ssurf[4];
1671
1672 // Select surface
1673 //
1674 G4double select = ssurf[5]*G4QuickRand();
1675 G4int k = 5;
1676 k -= (G4int)(select <= ssurf[4]);
1677 k -= (G4int)(select <= ssurf[3]);
1678 k -= (G4int)(select <= ssurf[2]);
1679 k -= (G4int)(select <= ssurf[1]);
1680 k -= (G4int)(select <= ssurf[0]);
1681
1682 // Generate point on selected surface
1683 //
1684 switch(k)
1685 {
1686 case 0: // start phi cut
1687 {
1688 G4double r = Rmin + (Rmax - Rmin)*G4Qui
1689 return { r*cosSPhi, r*sinSPhi, hz*G4Qui
1690 }
1691 case 1: // end phi cut
1692 {
1693 G4double r = Rmin + (Rmax - Rmin)*G4Qui
1694 return { r*cosEPhi, r*sinEPhi, hz*G4Qui
1695 }
1696 case 2: // base at -dz
1697 {
1698 G4double r = std::sqrt(RRmin + (RRmax -
1699 G4double phi = fSPhi + fDPhi*G4QuickRan
1700 return { r*std::cos(phi), r*std::sin(ph
1701 }
1702 case 3: // base at +dz
1703 {
1704 G4double r = std::sqrt(RRmin + (RRmax -
1705 G4double phi = fSPhi + fDPhi*G4QuickRan
1706 return { r*std::cos(phi), r*std::sin(ph
1707 }
1708 case 4: // external lateral surface
1709 {
1710 G4double phi = fSPhi + fDPhi*G4QuickRan
1711 G4double z = hz*G4QuickRand() - fDz;
1712 G4double x = Rmax*std::cos(phi);
1713 G4double y = Rmax*std::sin(phi);
1714 return { x,y,z };
1715 }
1716 case 5: // internal lateral surface
1717 {
1718 G4double phi = fSPhi + fDPhi*G4QuickRan
1719 G4double z = hz*G4QuickRand() - fDz;
1720 G4double x = Rmin*std::cos(phi);
1721 G4double y = Rmin*std::sin(phi);
1722 return { x,y,z };
1723 }
1724 }
1725 return {0., 0., 0.};
1726 }
1727
1728 /////////////////////////////////////////////
1729 //
1730 // Methods for visualisation
1731
1732 void G4Tubs::DescribeYourselfTo ( G4VGraphics
1733 {
1734 scene.AddSolid (*this) ;
1735 }
1736
1737 G4Polyhedron* G4Tubs::CreatePolyhedron () con
1738 {
1739 return new G4PolyhedronTubs (fRMin, fRMax,
1740 }
1741
1742 #endif
1743