| Geant4 Cross Reference |
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 // G4Tubs implementation
27 //
28 // 1994-95 P.Kent: first implementation
29 // 08.08.00 V.Grichine: more stable roots of 2-equation in DistanceToOut(p,v,..)
30 // 07.12.00 V.Grichine: phi-section algorithm was changed in Inside(p)
31 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J.Apostolakis proposal
32 // 24.08.16 E.Tcherniaev: reimplemented CalculateExtent().
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 angles so 0<sphi+dpshi<=2_PI
56 // - note if pdphi>2PI then reset to 2PI
57
58 G4Tubs::G4Tubs( const G4String& pName,
59 G4double pRMin, G4double pRMax,
60 G4double pDz,
61 G4double pSPhi, G4double pDPhi )
62 : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz),
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::GetInstance()->GetRadialTolerance();
68 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
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 (" << pDz << ") in solid: " << GetName();
78 G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
79 }
80 if ( (pRMin >= pRMax) || (pRMin < 0) ) // Check radii
81 {
82 std::ostringstream message;
83 message << "Invalid values for radii in solid: " << GetName()
84 << G4endl
85 << " pRMin = " << pRMin << ", pRMax = " << pRMax;
86 G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
87 }
88
89 // Check angles
90 //
91 CheckPhiAngles(pSPhi, pDPhi);
92 }
93
94 ///////////////////////////////////////////////////////////////////////
95 //
96 // Fake default constructor - sets only member data and allocates memory
97 // for usage restricted to object persistency.
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; kAngTolerance = rhs.kAngTolerance;
133 fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = rhs.fDz;
134 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
135 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; cosHDPhi = rhs.cosHDPhi;
136 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
137 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
138 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
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 replication mechanism dimension
152 // computation & modification.
153
154 void G4Tubs::ComputeDimensions( G4VPVParameterisation* p,
155 const G4int n,
156 const G4VPhysicalVolume* pRep )
157 {
158 p->ComputeDimensions(*this,n,pRep) ;
159 }
160
161 /////////////////////////////////////////////////////////////////////////
162 //
163 // Get bounding box
164
165 void G4Tubs::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
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(),GetCosStartPhi(),
178 GetSinEndPhi(),GetCosEndPhi(),
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.y() || pMin.z() >= pMax.z())
192 {
193 std::ostringstream message;
194 message << "Bad bounding box (min >= max) for solid: "
195 << GetName() << " !"
196 << "\npMin = " << pMin
197 << "\npMax = " << pMax;
198 G4Exception("G4Tubs::BoundingLimits()", "GeomMgt0001",
199 JustWarning, message);
200 DumpInfo();
201 }
202 }
203
204 /////////////////////////////////////////////////////////////////////////
205 //
206 // Calculate extent under transform and specified limit
207
208 G4bool G4Tubs::CalculateExtent( const EAxis pAxis,
209 const G4VoxelLimits& pVoxelLimit,
210 const G4AffineTransform& pTransform,
211 G4double& pMin,
212 G4double& pMax ) const
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,pVoxelLimit,pTransform,pMin,pMax);
224 #endif
225 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
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 extent
237 //
238 const G4int NSTEPS = 24; // number of steps for whole circle
239 G4double astep = twopi/NSTEPS; // max angle for one step
240 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
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 consists of two polygons,
250 // in other cases it is a sequence of quadrilaterals
251 if (rmin == 0 && dphi == twopi)
252 {
253 G4double sinCur = sinHalf;
254 G4double cosCur = cosHalf;
255
256 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
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 *> polygons(2);
267 polygons[0] = &baseA;
268 polygons[1] = &baseB;
269 G4BoundingEnvelope benv(bmin,bmax,polygons);
270 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
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 + cosStart*sinHalf;
279 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
280
281 // set quadrilaterals
282 G4ThreeVectorList pols[NSTEPS+2];
283 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
284 pols[0][0].set(rmin*cosStart,rmin*sinStart, dz);
285 pols[0][1].set(rmin*cosStart,rmin*sinStart,-dz);
286 pols[0][2].set(rmax*cosStart,rmax*sinStart,-dz);
287 pols[0][3].set(rmax*cosStart,rmax*sinStart, dz);
288 for (G4int k=1; k<ksteps+1; ++k)
289 {
290 pols[k][0].set(rmin*cosCur,rmin*sinCur, dz);
291 pols[k][1].set(rmin*cosCur,rmin*sinCur,-dz);
292 pols[k][2].set(rext*cosCur,rext*sinCur,-dz);
293 pols[k][3].set(rext*cosCur,rext*sinCur, dz);
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*sinEnd, dz);
300 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,-dz);
301 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,-dz);
302 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd, dz);
303
304 // set envelope and calculate extent
305 std::vector<const G4ThreeVectorList *> polygons;
306 polygons.resize(ksteps+2);
307 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
308 G4BoundingEnvelope benv(bmin,bmax,polygons);
309 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
310 }
311 return exist;
312 }
313
314 ///////////////////////////////////////////////////////////////////////////
315 //
316 // Return whether point inside/outside/on surface
317
318 EInside G4Tubs::Inside( const G4ThreeVector& p ) const
319 {
320 G4double r2,pPhi,tolRMin,tolRMax;
321 EInside in = kOutside ;
322
323 if (std::fabs(p.z()) <= fDz - halfCarTolerance)
324 {
325 r2 = p.x()*p.x() + p.y()*p.y() ;
326
327 if (fRMin != 0.0) { tolRMin = fRMin + halfRadTolerance ; }
328 else { tolRMin = 0 ; }
329
330 tolRMax = fRMax - halfRadTolerance ;
331
332 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
333 {
334 if ( fPhiFullTube )
335 {
336 in = kInside ;
337 }
338 else
339 {
340 // Try inner tolerant phi boundaries (=>inside)
341 // if not inside, try outer tolerant phi boundaries
342
343 if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
344 && (std::fabs(p.y())<=halfCarTolerance) )
345 {
346 in=kSurface;
347 }
348 else
349 {
350 pPhi = std::atan2(p.y(),p.x()) ;
351 if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi
352
353 if ( fSPhi >= 0 )
354 {
355 if ( (std::fabs(pPhi) < halfAngTolerance)
356 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
357 {
358 pPhi += twopi ; // 0 <= pPhi < 2pi
359 }
360 if ( (pPhi >= fSPhi + halfAngTolerance)
361 && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
362 {
363 in = kInside ;
364 }
365 else if ( (pPhi >= fSPhi - halfAngTolerance)
366 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
367 {
368 in = kSurface ;
369 }
370 }
371 else // fSPhi < 0
372 {
373 if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
374 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} //kOutside
375 else if ( (pPhi <= fSPhi + twopi + halfAngTolerance)
376 && (pPhi >= fSPhi + fDPhi - halfAngTolerance) )
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 <= tolRMax*tolRMax) )
396 {
397 if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
398 { // Continuous in phi or on z-axis
399 in = kSurface ;
400 }
401 else // Try outer tolerant phi boundaries only
402 {
403 pPhi = std::atan2(p.y(),p.x()) ;
404
405 if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi
406 if ( fSPhi >= 0 )
407 {
408 if ( (std::fabs(pPhi) < halfAngTolerance)
409 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
410 {
411 pPhi += twopi ; // 0 <= pPhi < 2pi
412 }
413 if ( (pPhi >= fSPhi - halfAngTolerance)
414 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
415 {
416 in = kSurface ;
417 }
418 }
419 else // fSPhi < 0
420 {
421 if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
422 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside
423 else
424 {
425 in = kSurface ;
426 }
427 }
428 }
429 }
430 }
431 }
432 else if (std::fabs(p.z()) <= fDz + halfCarTolerance)
433 { // Check within tolerant r limits
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 <= tolRMax*tolRMax) )
441 {
442 if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
443 { // Continuous in phi or on z-axis
444 in = kSurface ;
445 }
446 else // Try outer tolerant phi boundaries
447 {
448 pPhi = std::atan2(p.y(),p.x()) ;
449
450 if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi
451 if ( fSPhi >= 0 )
452 {
453 if ( (std::fabs(pPhi) < halfAngTolerance)
454 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
455 {
456 pPhi += twopi ; // 0 <= pPhi < 2pi
457 }
458 if ( (pPhi >= fSPhi - halfAngTolerance)
459 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
460 {
461 in = kSurface;
462 }
463 }
464 else // fSPhi < 0
465 {
466 if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
467 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
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 divided sides
483 // - unsafe if point close to z axis a rmin=0 - no explicit checks
484
485 G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p ) const
486 {
487 G4int noSurfaces = 0;
488 G4double rho, pPhi;
489 G4double distZ, distRMin, distRMax;
490 G4double distSPhi = kInfinity, distEPhi = kInfinity;
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 (0,0,z)
503 {
504 if ( rho > halfCarTolerance )
505 {
506 pPhi = std::atan2(p.y(),p.x());
507
508 if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; }
509 else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; }
510
511 distSPhi = std::fabs( pPhi - fSPhi );
512 distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
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 = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
523
524 if( distRMax <= halfCarTolerance )
525 {
526 ++noSurfaces;
527 sumnorm += nR;
528 }
529 if( (fRMin != 0.0) && (distRMin <= halfCarTolerance) )
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)", "GeomSolids1002",
557 JustWarning, "Point p is not on surface !?" );
558 G4long oldprc = G4cout.precision(20);
559 G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
560 << G4endl << G4endl;
561 G4cout.precision(oldprc) ;
562 #endif
563 norm = ApproxSurfaceNormal(p);
564 }
565 else if ( noSurfaces == 1 ) { norm = sumnorm; }
566 else { norm = sumnorm.unit(); }
567
568 return norm;
569 }
570
571 /////////////////////////////////////////////////////////////////////////////
572 //
573 // Algorithm for SurfaceNormal() following the original specification
574 // for points not on the surface
575
576 G4ThreeVector G4Tubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const
577 {
578 ENorm side ;
579 G4ThreeVector norm ;
580 G4double rho, phi ;
581 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
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) ) // Protected against (0,0,z)
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 + twopi))*rho ;
624 }
625 else
626 {
627 distSPhi = std::fabs(phi - fSPhi)*rho ;
628 }
629 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
630
631 if (distSPhi < distEPhi) // Find new minimum
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()/rho, 0) ;
651 break ;
652 }
653 case kNRMax : // Outer radius
654 {
655 norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
656 break ;
657 }
658 case kNZ : // + or - dz
659 {
660 if ( p.z() > 0 ) { norm = G4ThreeVector(0,0,1) ; }
661 else { norm = G4ThreeVector(0,0,-1); }
662 break ;
663 }
664 case kNSPhi:
665 {
666 norm = G4ThreeVector(sinSPhi, -cosSPhi, 0) ;
667 break ;
668 }
669 case kNEPhi:
670 {
671 norm = G4ThreeVector(-sinEPhi, cosEPhi, 0) ;
672 break;
673 }
674 default: // Should never reach this case ...
675 {
676 DumpInfo();
677 G4Exception("G4Tubs::ApproxSurfaceNormal()",
678 "GeomSolids1002", JustWarning,
679 "Undefined side for valid surface normal to solid.");
680 break ;
681 }
682 }
683 return norm;
684 }
685
686 ////////////////////////////////////////////////////////////////////
687 //
688 //
689 // Calculate distance to shape from outside, along normalised vector
690 // - return kInfinity if no intersection, or intersection distance <= tolerance
691 //
692 // - Compute the intersection with the z planes
693 // - if at valid r, phi, return
694 //
695 // -> If point is outer outer radius, compute intersection with rmax
696 // - if at valid phi,z return
697 //
698 // -> Compute intersection with inner radius, taking largest +ve root
699 // - if valid (in z,phi), save intersction
700 //
701 // -> If phi segmented, compute intersections with phi half planes
702 // - return smallest of valid phi intersections and
703 // inner radius intersection
704 //
705 // NOTE:
706 // - 'if valid' implies tolerant checking of intersection points
707
708 G4double G4Tubs::DistanceToIn( const G4ThreeVector& p,
709 const G4ThreeVector& v ) const
710 {
711 G4double snxt = kInfinity ; // snxt = default return value
712 G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared
713 G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
714 const G4double dRmax = 100.*fRMax;
715
716 // Intersection point variables
717 //
718 G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
719 G4double t1, t2, t3, b, c, d ; // Quadratic solver variables
720
721 // Calculate tolerant rmin and rmax
722
723 if (fRMin > kRadTolerance)
724 {
725 tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
726 tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
727 }
728 else
729 {
730 tolORMin2 = 0.0 ;
731 tolIRMin2 = 0.0 ;
732 }
733 tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
734 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
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 in -Z or visa versa
744 {
745 sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance
746
747 if(sd < 0.0) { sd = 0.0; }
748
749 xi = p.x() + sd*v.x() ; // Intersection coords
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 <= tolIRMax2))
756 {
757 if (!fPhiFullTube && (rho2 != 0.0))
758 {
759 // Psi = angle made with central (average) phi of shape
760 //
761 inum = xi*cosCPhi + yi*sinCPhi ;
762 iden = std::sqrt(rho2) ;
763 cosPsi = inum/iden ;
764 if (cosPsi >= cosHDPhiIT) { return sd ; }
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 heading away
776 // -> cannot intersect
777 }
778 }
779
780 // -> Can not intersect z surfaces
781 //
782 // Intersection with rmax (possible return) and rmin (must also check phi)
783 //
784 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
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.y)+p.x^2+p.y^2-R^2=0
789 // t1 t2 t3
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 axis
796 {
797 b = t2/t1 ;
798 c = t3 - fRMax*fRMax ;
799 if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case
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 errors due to precision issues on
813 { // 64 bits systems. Split long distances and recompute
814 G4double fTerm = sd-std::fmod(sd,dRmax);
815 sd = fTerm + DistanceToIn(p+fTerm*v,v);
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 reqd
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*sinCPhi)/fRMax ;
833 if (cosPsi >= cosHDPhiIT) { return sd ; }
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 tubs (-> 0 to in)
843
844 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
845 {
846 // Inside both radii, delta r -ve, inside z extent
847
848 if (!fPhiFullTube)
849 {
850 inum = p.x()*cosCPhi + p.y()*sinCPhi ;
851 iden = std::sqrt(t3) ;
852 cosPsi = inum/iden ;
853 if (cosPsi >= cosHDPhiIT)
854 {
855 // In the old version, the small negative tangent for the point
856 // on surface was not taken in account, and returning 0.0 ...
857 // New version: check the tangent for the point on surface and
858 // if no intersection, return kInfinity, if intersection instead
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)); // using safe solution
873 // for quadratic equation
874 if ( snxt < halfCarTolerance ) { snxt=0; }
875 return snxt ;
876 }
877 else
878 {
879 return kInfinity;
880 }
881 }
882 }
883 }
884 else
885 {
886 // In the old version, the small negative tangent for the point
887 // on surface was not taken in account, and returning 0.0 ...
888 // New version: check the tangent for the point on surface and
889 // if no intersection, return kInfinity, if intersection instead
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)); // using safe solution
904 // for quadratic equation
905 if ( snxt < halfCarTolerance ) { snxt=0; }
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 cylinder intersection
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 outside and know rmax Hit was bad
923 // - If on surface of rmin also need farthest root
924
925 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
926 if (sd >= -halfCarTolerance) // check forwards
927 {
928 // Check z intersection
929 //
930 if(sd < 0.0) { sd = 0.0; }
931 if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen
932 { // 64 bits systems. Split long distances and recompute
933 G4double fTerm = sd-std::fmod(sd,dRmax);
934 sd = fTerm + DistanceToIn(p+fTerm*v,v);
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*sinCPhi)*fInvRmin;
950 if (cosPsi >= cosHDPhiIT)
951 {
952 // Good inner radius isect
953 // - but earlier phi isect still possible
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 up to kCarTolerance*0.5
967 //
968 // o NOTE: Large duplication of code between sphi & ephi checks
969 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
970 // intersection check <=0 -> >=0
971 // -> use some form of loop Construct ?
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 normal dirn
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) && (rho2 <= tolIRMax2) )
998 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
999 && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 )
1000 && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) )
1001 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1002 && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1003 && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) )
1004 {
1005 // z and r intersections good
1006 // - check intersecting with correct half-plane
1007 //
1008 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
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 normal dirn
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) && (rho2 <= tolIRMax2) )
1037 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1038 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1039 && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1040 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1041 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1042 && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1043 {
1044 // z and r intersections good
1045 // - check intersecting with correct half-plane
1046 //
1047 if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = sd; }
1048 } //?? >=-halfCarTolerance
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, along normalised vector
1061 // - return kInfinity if no intersection, or intersection distance <= tolerance
1062 //
1063 // - Compute the intersection with the z planes
1064 // - if at valid r, phi, return
1065 //
1066 // -> If point is outer outer radius, compute intersection with rmax
1067 // - if at valid phi,z return
1068 //
1069 // -> Compute intersection with inner radius, taking largest +ve root
1070 // - if valid (in z,phi), save intersction
1071 //
1072 // -> If phi segmented, compute intersections with phi half planes
1073 // - return smallest of valid phi intersections and
1074 // inner radius intersection
1075 //
1076 // NOTE:
1077 // - Precalculations for phi trigonometry are Done `just in time'
1078 // - `if valid' implies tolerant checking of intersection points
1079 // Calculate distance (<= actual) to closest surface of shape from outside
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 G4ThreeVector& p ) const
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)/rho ;
1103
1104 if ( cosPsi < cosHDPhi )
1105 {
1106 // Point lies outside phi range
1107
1108 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1109 {
1110 safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1111 }
1112 else
1113 {
1114 safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
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 from `inside', allowing for tolerance
1126 // - Only Calc rmax intersection if no valid rmin intersection
1127
1128 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p,
1129 const G4ThreeVector& v,
1130 const G4bool calcNorm,
1131 G4bool* validNorm,
1132 G4ThreeVector* n ) const
1133 {
1134 ESide side=kNull , sider=kNull, sidephi=kNull ;
1135 G4double snxt, srd=kInfinity, sphi=kInfinity, pdist ;
1136 G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1137
1138 // Vars for phi intersection:
1139
1140 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
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 perpendicular to z axis
1184 side = kNull;
1185 }
1186
1187 // Radial Intersections
1188 //
1189 // Find intersection with cylinders at rmax/rmin
1190 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
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.y)+p.x^2+p.y^2-R^2=0
1195 //
1196 // t1 t2 t3
1197
1198 t1 = 1.0 - v.z()*v.z() ; // since v normalised
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*fRMax*fRMax; }
1203 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz
1204
1205 if ( t1 > 0 ) // Check not parallel
1206 {
1207 // Calculate srd, r exit distance
1208
1209 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1210 {
1211 // Delta r not negative => leaving via rmax
1212
1213 deltaR = t3 - fRMax*fRMax ;
1214
1215 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
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::sqrt(d2)); }
1224 else { srd = 0.; }
1225 sider = kRMax ;
1226 }
1227 else
1228 {
1229 // On tolerant boundary & heading outwards (or perpendicular to)
1230 // outer radial surface -> leaving immediately
1231
1232 if ( calcNorm )
1233 {
1234 G4double invRho = FastInverseRxy( p, fInvRmax, kNormTolerance );
1235 *n = G4ThreeVector(p.x()*invRho,p.y()*invRho,0) ;
1236 *validNorm = true ;
1237 }
1238 return snxt = 0 ; // Leaving by rmax immediately
1239 }
1240 }
1241 else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection
1242 {
1243 roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1244
1245 if ( (fRMin != 0.0) && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
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>kRadTolerance*0.5
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 be rmax intersect
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<kRadTolerance
1281 // (v is perpendicular to the surface)
1282 {
1283 if (calcNorm)
1284 {
1285 G4double invRho = FastInverseRxy( p, fInvRmax, kNormTolerance );
1286 *n = G4ThreeVector(p.x()*invRho,p.y()*invRho,0) ;
1287 *validNorm = true ;
1288 }
1289 return snxt = 0.0;
1290 }
1291 }
1292 }
1293 else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1294 // No rmin intersect -> must be rmax intersect
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<kRadTolerance
1306 // (v is perpendicular to the surface)
1307 {
1308 if (calcNorm)
1309 {
1310 G4double invRho = FastInverseRxy( p, fInvRmax, kNormTolerance );
1311 *n = G4ThreeVector(p.x()*invRho,p.y()*invRho,0) ;
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 correction
1324 // of the difference in domain of atan2 and Sphi
1325 //
1326 vphi = std::atan2(v.y(),v.x()) ;
1327
1328 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1329 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1330
1331
1332 if ( (p.x() != 0.0) || (p.y() != 0.0) ) // Check if on z axis (rho not needed later)
1333 {
1334 // pDist -ve when inside
1335
1336 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1337 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1338
1339 // Comp -ve when in direction of outwards normal
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 <= halfCarTolerance)
1347 && (pDistE <= halfCarTolerance) ) )
1348 || ( (fDPhi > pi) && ((pDistS <= halfCarTolerance)
1349 || (pDistE <= halfCarTolerance) ) ) )
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 correct half-plane
1363 // (if not -> no intersect)
1364 //
1365 if((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1366 {
1367 sidephi = kSPhi;
1368 if (((fSPhi-halfAngTolerance)<=vphi)
1369 &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1370 {
1371 sphi = kInfinity;
1372 }
1373 }
1374 else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1375 {
1376 sphi = kInfinity ;
1377 }
1378 else
1379 {
1380 sidephi = kSPhi ;
1381 if ( pDistS > -halfCarTolerance )
1382 {
1383 sphi = 0.0 ; // Leave by sphi immediately
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 < starting phi intersection
1402 //
1403 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1404 {
1405 xi = p.x() + sphi2*v.x() ;
1406 yi = p.y() + sphi2*v.y() ;
1407
1408 if((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1409 {
1410 // Leaving via ending phi
1411 //
1412 if( (fSPhi-halfAngTolerance > vphi)
1413 ||(fSPhi+fDPhi+halfAngTolerance < vphi) )
1414 {
1415 sidephi = kEPhi ;
1416 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1417 else { sphi = 0.0 ; }
1418 }
1419 }
1420 else // Check intersecting with correct half-plane
1421
1422 if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1423 {
1424 // Leaving via ending phi
1425 //
1426 sidephi = kEPhi ;
1427 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1428 else { sphi = 0.0 ; }
1429 }
1430 }
1431 }
1432 }
1433 else
1434 {
1435 sphi = kInfinity ;
1436 }
1437 }
1438 else
1439 {
1440 // On z axis + travel not || to z axis -> if phi of vector direction
1441 // within phi of shape, Step limited by rmax, else Step =0
1442
1443 if ( (fSPhi - halfAngTolerance <= vphi)
1444 && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1445 {
1446 sphi = kInfinity ;
1447 }
1448 else
1449 {
1450 sidephi = kSPhi ; // arbitrary
1451 sphi = 0.0 ;
1452 }
1453 }
1454 if (sphi < snxt) // Order intersecttions
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 normalised
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,0) ;
1477 *validNorm = true ;
1478 break ;
1479
1480 case kRMin:
1481 *validNorm = false ; // Rmin is inconvex
1482 break ;
1483
1484 case kSPhi:
1485 if ( fDPhi <= pi )
1486 {
1487 *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
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,0) ;
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 surface normal to solid."
1524 << G4endl
1525 << "Position:" << G4endl << G4endl
1526 << "p.x() = " << p.x()/mm << " mm" << G4endl
1527 << "p.y() = " << p.y()/mm << " mm" << G4endl
1528 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1529 << "Direction:" << G4endl << G4endl
1530 << "v.x() = " << v.x() << G4endl
1531 << "v.y() = " << v.y() << G4endl
1532 << "v.z() = " << v.z() << G4endl << G4endl
1533 << "Proposed distance :" << G4endl << G4endl
1534 << "snxt = " << snxt/mm << " mm" << G4endl ;
1535 message.precision(oldprc) ;
1536 G4Exception("G4Tubs::DistanceToOut(p,v,..)", "GeomSolids1002",
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 surface of shape from inside
1549
1550 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const
1551 {
1552 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
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 << G4endl ;
1562 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1563 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1564 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1565 G4cout.precision(oldprc) ;
1566 G4Exception("G4Tubs::DistanceToOut(p)", "GeomSolids1002",
1567 JustWarning, "Point p is outside !?");
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 closest phi plane
1588 //
1589 if ( !fPhiFullTube )
1590 {
1591 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1592 {
1593 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1594 }
1595 else
1596 {
1597 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
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::ostream& os ) const
1629 {
1630 G4long oldprc = os.precision(16);
1631 os << "-----------------------------------------------------------\n"
1632 << " *** Dump for solid - " << GetName() << " ***\n"
1633 << " ===================================================\n"
1634 << " Solid type: G4Tubs\n"
1635 << " Parameters: \n"
1636 << " inner radius : " << fRMin/mm << " mm \n"
1637 << " outer radius : " << fRMax/mm << " mm \n"
1638 << " half length Z: " << fDz/mm << " mm \n"
1639 << " starting phi : " << fSPhi/degree << " degrees \n"
1640 << " delta phi : " << fDPhi/degree << " degrees \n"
1641 << "-----------------------------------------------------------\n";
1642 os.precision(oldprc);
1643
1644 return os;
1645 }
1646
1647 /////////////////////////////////////////////////////////////////////////
1648 //
1649 // GetPointOnSurface
1650
1651 G4ThreeVector G4Tubs::GetPointOnSurface() const
1652 {
1653 G4double Rmax = fRMax;
1654 G4double Rmin = fRMin;
1655 G4double hz = 2.*fDz; // height
1656 G4double lext = fDPhi*Rmax; // length of external circular arc
1657 G4double lint = fDPhi*Rmin; // length of internal circular arc
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*(Rmax - Rmin);
1665 G4double ssurf[6] = { scut, scut, sbase, sbase, hz*lext, hz*lint };
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)*G4QuickRand();
1689 return { r*cosSPhi, r*sinSPhi, hz*G4QuickRand() - fDz };
1690 }
1691 case 1: // end phi cut
1692 {
1693 G4double r = Rmin + (Rmax - Rmin)*G4QuickRand();
1694 return { r*cosEPhi, r*sinEPhi, hz*G4QuickRand() - fDz };
1695 }
1696 case 2: // base at -dz
1697 {
1698 G4double r = std::sqrt(RRmin + (RRmax - RRmin)*G4QuickRand());
1699 G4double phi = fSPhi + fDPhi*G4QuickRand();
1700 return { r*std::cos(phi), r*std::sin(phi), -fDz };
1701 }
1702 case 3: // base at +dz
1703 {
1704 G4double r = std::sqrt(RRmin + (RRmax - RRmin)*G4QuickRand());
1705 G4double phi = fSPhi + fDPhi*G4QuickRand();
1706 return { r*std::cos(phi), r*std::sin(phi), fDz };
1707 }
1708 case 4: // external lateral surface
1709 {
1710 G4double phi = fSPhi + fDPhi*G4QuickRand();
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*G4QuickRand();
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 ( G4VGraphicsScene& scene ) const
1733 {
1734 scene.AddSolid (*this) ;
1735 }
1736
1737 G4Polyhedron* G4Tubs::CreatePolyhedron () const
1738 {
1739 return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi) ;
1740 }
1741
1742 #endif
1743