Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4Cons.cc
Go to the documentation of this file.
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$
28// GEANT4 tag $Name: $
29//
30//
31// class G4Cons
32//
33// Implementation for G4Cons class
34//
35// History:
36//
37// 05.04.12 M.Kelsey: GetPointOnSurface() throw flat in sqrt(r)
38// 12.10.09 T.Nikitina: Added to DistanceToIn(p,v) check on the direction in
39// case of point on surface
40// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
41// 13.09.96 V.Grichine: Review and final modifications
42// ~1994 P.Kent: Created, as main part of the geometry prototype
43// --------------------------------------------------------------------
44
45#include "G4Cons.hh"
46
47#include "G4VoxelLimits.hh"
48#include "G4AffineTransform.hh"
49#include "G4GeometryTolerance.hh"
50
51#include "G4VPVParameterisation.hh"
52
53#include "meshdefs.hh"
54
55#include "Randomize.hh"
56
57#include "G4VGraphicsScene.hh"
58#include "G4Polyhedron.hh"
59#include "G4NURBS.hh"
60#include "G4NURBSbox.hh"
61
62using namespace CLHEP;
63
64////////////////////////////////////////////////////////////////////////
65//
66// Private enum: Not for external use - used by distanceToOut
67
68enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
69
70// used by normal
71
72enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
73
74//////////////////////////////////////////////////////////////////////////
75//
76// constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
77// - note if pDPhi>2PI then reset to 2PI
78
79G4Cons::G4Cons( const G4String& pName,
80 G4double pRmin1, G4double pRmax1,
81 G4double pRmin2, G4double pRmax2,
82 G4double pDz,
83 G4double pSPhi, G4double pDPhi)
84 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
85 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
86{
87 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
88 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
89
90 // Check z-len
91 //
92 if ( pDz < 0 )
93 {
94 std::ostringstream message;
95 message << "Invalid Z half-length for Solid: " << GetName() << G4endl
96 << " hZ = " << pDz;
97 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
98 FatalException, message);
99 }
100
101 // Check radii
102 //
103 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
104 {
105 std::ostringstream message;
106 message << "Invalid values of radii for Solid: " << GetName() << G4endl
107 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
108 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
109 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
110 FatalException, message) ;
111 }
112 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
113 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
114
115 // Check angles
116 //
117 CheckPhiAngles(pSPhi, pDPhi);
118}
119
120///////////////////////////////////////////////////////////////////////
121//
122// Fake default constructor - sets only member data and allocates memory
123// for usage restricted to object persistency.
124//
125G4Cons::G4Cons( __void__& a )
126 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
127 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
128 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
129 cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
130 fPhiFullCone(false)
131{
132}
133
134///////////////////////////////////////////////////////////////////////
135//
136// Destructor
137
138G4Cons::~G4Cons()
139{
140}
141
142//////////////////////////////////////////////////////////////////////////
143//
144// Copy constructor
145
146G4Cons::G4Cons(const G4Cons& rhs)
147 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
148 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
149 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
150 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
151 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
152 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
153 cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone)
154{
155}
156
157//////////////////////////////////////////////////////////////////////////
158//
159// Assignment operator
160
161G4Cons& G4Cons::operator = (const G4Cons& rhs)
162{
163 // Check assignment to self
164 //
165 if (this == &rhs) { return *this; }
166
167 // Copy base class data
168 //
169 G4CSGSolid::operator=(rhs);
170
171 // Copy data
172 //
173 kRadTolerance = rhs.kRadTolerance;
174 kAngTolerance = rhs.kAngTolerance;
175 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
176 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
177 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
178 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
179 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
180 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
181 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
182 fPhiFullCone = rhs.fPhiFullCone;
183
184 return *this;
185}
186
187/////////////////////////////////////////////////////////////////////
188//
189// Return whether point inside/outside/on surface
190
191EInside G4Cons::Inside(const G4ThreeVector& p) const
192{
193 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
194 EInside in;
195 static const G4double halfCarTolerance=kCarTolerance*0.5;
196 static const G4double halfRadTolerance=kRadTolerance*0.5;
197 static const G4double halfAngTolerance=kAngTolerance*0.5;
198
199 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
200 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
201 else { in = kInside; }
202
203 r2 = p.x()*p.x() + p.y()*p.y() ;
204 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
205 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
206
207 // rh2 = rh*rh;
208
209 tolRMin = rl - halfRadTolerance;
210 if ( tolRMin < 0 ) { tolRMin = 0; }
211 tolRMax = rh + halfRadTolerance;
212
213 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
214
215 if (rl) { tolRMin = rl + halfRadTolerance; }
216 else { tolRMin = 0.0; }
217 tolRMax = rh - halfRadTolerance;
218
219 if (in == kInside) // else it's kSurface already
220 {
221 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
222 }
223 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
224 {
225 pPhi = std::atan2(p.y(),p.x()) ;
226
227 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
228 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
229
230 if ( (pPhi < fSPhi - halfAngTolerance) ||
231 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
232
233 else if (in == kInside) // else it's kSurface anyway already
234 {
235 if ( (pPhi < fSPhi + halfAngTolerance) ||
236 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
237 }
238 }
239 else if ( !fPhiFullCone ) { in = kSurface; }
240
241 return in ;
242}
243
244/////////////////////////////////////////////////////////////////////////
245//
246// Dispatch to parameterisation for replication mechanism dimension
247// computation & modification.
248
249void G4Cons::ComputeDimensions( G4VPVParameterisation* p,
250 const G4int n,
251 const G4VPhysicalVolume* pRep )
252{
253 p->ComputeDimensions(*this,n,pRep) ;
254}
255
256
257///////////////////////////////////////////////////////////////////////////
258//
259// Calculate extent under transform and specified limit
260
261G4bool G4Cons::CalculateExtent( const EAxis pAxis,
262 const G4VoxelLimits& pVoxelLimit,
263 const G4AffineTransform& pTransform,
264 G4double& pMin,
265 G4double& pMax ) const
266{
267 if ( !pTransform.IsRotated() && (fDPhi == twopi)
268 && (fRmin1 == 0) && (fRmin2 == 0) )
269 {
270 // Special case handling for unrotated solid cones
271 // Compute z/x/y mins and maxs for bounding box respecting limits,
272 // with early returns if outside limits. Then switch() on pAxis,
273 // and compute exact x and y limit for x/y case
274
275 G4double xoffset, xMin, xMax ;
276 G4double yoffset, yMin, yMax ;
277 G4double zoffset, zMin, zMax ;
278
279 G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
280 G4double xoff1, xoff2, yoff1, yoff2 ;
281
282 zoffset = pTransform.NetTranslation().z();
283 zMin = zoffset - fDz ;
284 zMax = zoffset + fDz ;
285
286 if (pVoxelLimit.IsZLimited())
287 {
288 if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) ||
289 (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
290 {
291 return false ;
292 }
293 else
294 {
295 if ( zMin < pVoxelLimit.GetMinZExtent() )
296 {
297 zMin = pVoxelLimit.GetMinZExtent() ;
298 }
299 if ( zMax > pVoxelLimit.GetMaxZExtent() )
300 {
301 zMax = pVoxelLimit.GetMaxZExtent() ;
302 }
303 }
304 }
305 xoffset = pTransform.NetTranslation().x() ;
306 RMax = (fRmax2 >= fRmax1) ? zMax : zMin ;
307 xMax = xoffset + (fRmax1 + fRmax2)*0.5 +
308 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
309 xMin = 2*xoffset-xMax ;
310
311 if (pVoxelLimit.IsXLimited())
312 {
313 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) ||
314 (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
315 {
316 return false ;
317 }
318 else
319 {
320 if ( xMin < pVoxelLimit.GetMinXExtent() )
321 {
322 xMin = pVoxelLimit.GetMinXExtent() ;
323 }
324 if ( xMax > pVoxelLimit.GetMaxXExtent() )
325 {
326 xMax=pVoxelLimit.GetMaxXExtent() ;
327 }
328 }
329 }
330 yoffset = pTransform.NetTranslation().y() ;
331 yMax = yoffset + (fRmax1 + fRmax2)*0.5 +
332 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
333 yMin = 2*yoffset-yMax ;
334 RMax = yMax - yoffset ; // = max radius due to Zmax/Zmin cuttings
335
336 if (pVoxelLimit.IsYLimited())
337 {
338 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) ||
339 (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
340 {
341 return false ;
342 }
343 else
344 {
345 if ( yMin < pVoxelLimit.GetMinYExtent() )
346 {
347 yMin = pVoxelLimit.GetMinYExtent() ;
348 }
349 if ( yMax > pVoxelLimit.GetMaxYExtent() )
350 {
351 yMax = pVoxelLimit.GetMaxYExtent() ;
352 }
353 }
354 }
355 switch (pAxis) // Known to cut cones
356 {
357 case kXAxis:
358 yoff1 = yoffset - yMin ;
359 yoff2 = yMax - yoffset ;
360
361 if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
362 { // => no change
363 pMin = xMin ;
364 pMax = xMax ;
365 }
366 else
367 {
368 // Y limits don't cross max/min x => compute max delta x,
369 // hence new mins/maxs
370 delta=RMax*RMax-yoff1*yoff1;
371 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
372 delta=RMax*RMax-yoff2*yoff2;
373 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
374 maxDiff = (diff1>diff2) ? diff1:diff2 ;
375 newMin = xoffset - maxDiff ;
376 newMax = xoffset + maxDiff ;
377 pMin = ( newMin < xMin ) ? xMin : newMin ;
378 pMax = ( newMax > xMax) ? xMax : newMax ;
379 }
380 break ;
381
382 case kYAxis:
383 xoff1 = xoffset - xMin ;
384 xoff2 = xMax - xoffset ;
385
386 if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
387 { // => no change
388 pMin = yMin ;
389 pMax = yMax ;
390 }
391 else
392 {
393 // X limits don't cross max/min y => compute max delta y,
394 // hence new mins/maxs
395 delta=RMax*RMax-xoff1*xoff1;
396 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
397 delta=RMax*RMax-xoff2*xoff2;
398 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
399 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
400 newMin = yoffset - maxDiff ;
401 newMax = yoffset + maxDiff ;
402 pMin = (newMin < yMin) ? yMin : newMin ;
403 pMax = (newMax > yMax) ? yMax : newMax ;
404 }
405 break ;
406
407 case kZAxis:
408 pMin = zMin ;
409 pMax = zMax ;
410 break ;
411
412 default:
413 break ;
414 }
415 pMin -= kCarTolerance ;
416 pMax += kCarTolerance ;
417
418 return true ;
419 }
420 else // Calculate rotated vertex coordinates
421 {
422 G4int i, noEntries, noBetweenSections4 ;
423 G4bool existsAfterClip = false ;
424 G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
425
426 pMin = +kInfinity ;
427 pMax = -kInfinity ;
428
429 noEntries = vertices->size() ;
430 noBetweenSections4 = noEntries-4 ;
431
432 for ( i = 0 ; i < noEntries ; i += 4 )
433 {
434 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
435 }
436 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
437 {
438 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
439 }
440 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
441 {
442 existsAfterClip = true ;
443
444 // Add 2*tolerance to avoid precision troubles
445
446 pMin -= kCarTolerance ;
447 pMax += kCarTolerance ;
448 }
449 else
450 {
451 // Check for case where completely enveloping clipping volume
452 // If point inside then we are confident that the solid completely
453 // envelopes the clipping volume. Hence set min/max extents according
454 // to clipping volume extents along the specified axis.
455
456 G4ThreeVector clipCentre(
457 (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
458 (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
459 (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5 ) ;
460
461 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
462 {
463 existsAfterClip = true ;
464 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
465 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
466 }
467 }
468 delete vertices ;
469 return existsAfterClip ;
470 }
471}
472
473////////////////////////////////////////////////////////////////////////
474//
475// Return unit normal of surface closest to p
476// - note if point on z axis, ignore phi divided sides
477// - unsafe if point close to z axis a rmin=0 - no explicit checks
478
479G4ThreeVector G4Cons::SurfaceNormal( const G4ThreeVector& p) const
480{
481 G4int noSurfaces = 0;
482 G4double rho, pPhi;
483 G4double distZ, distRMin, distRMax;
484 G4double distSPhi = kInfinity, distEPhi = kInfinity;
485 G4double tanRMin, secRMin, pRMin, widRMin;
486 G4double tanRMax, secRMax, pRMax, widRMax;
487
488 static const G4double delta = 0.5*kCarTolerance;
489 static const G4double dAngle = 0.5*kAngTolerance;
490
491 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
492 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
493
494 distZ = std::fabs(std::fabs(p.z()) - fDz);
495 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
496
497 tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
498 secRMin = std::sqrt(1 + tanRMin*tanRMin);
499 pRMin = rho - p.z()*tanRMin;
500 widRMin = fRmin2 - fDz*tanRMin;
501 distRMin = std::fabs(pRMin - widRMin)/secRMin;
502
503 tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
504 secRMax = std::sqrt(1+tanRMax*tanRMax);
505 pRMax = rho - p.z()*tanRMax;
506 widRMax = fRmax2 - fDz*tanRMax;
507 distRMax = std::fabs(pRMax - widRMax)/secRMax;
508
509 if (!fPhiFullCone) // Protected against (0,0,z)
510 {
511 if ( rho )
512 {
513 pPhi = std::atan2(p.y(),p.x());
514
515 if (pPhi < fSPhi-delta) { pPhi += twopi; }
516 else if (pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
517
518 distSPhi = std::fabs( pPhi - fSPhi );
519 distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
520 }
521 else if( !(fRmin1) || !(fRmin2) )
522 {
523 distSPhi = 0.;
524 distEPhi = 0.;
525 }
526 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
527 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
528 }
529 if ( rho > delta )
530 {
531 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
532 if (fRmin1 || fRmin2)
533 {
534 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
535 }
536 }
537
538 if( distRMax <= delta )
539 {
540 noSurfaces ++;
541 sumnorm += nR;
542 }
543 if( (fRmin1 || fRmin2) && (distRMin <= delta) )
544 {
545 noSurfaces ++;
546 sumnorm += nr;
547 }
548 if( !fPhiFullCone )
549 {
550 if (distSPhi <= dAngle)
551 {
552 noSurfaces ++;
553 sumnorm += nPs;
554 }
555 if (distEPhi <= dAngle)
556 {
557 noSurfaces ++;
558 sumnorm += nPe;
559 }
560 }
561 if (distZ <= delta)
562 {
563 noSurfaces ++;
564 if ( p.z() >= 0.) { sumnorm += nZ; }
565 else { sumnorm -= nZ; }
566 }
567 if ( noSurfaces == 0 )
568 {
569#ifdef G4CSGDEBUG
570 G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
571 JustWarning, "Point p is not on surface !?" );
572#endif
573 norm = ApproxSurfaceNormal(p);
574 }
575 else if ( noSurfaces == 1 ) { norm = sumnorm; }
576 else { norm = sumnorm.unit(); }
577
578 return norm ;
579}
580
581////////////////////////////////////////////////////////////////////////////
582//
583// Algorithm for SurfaceNormal() following the original specification
584// for points not on the surface
585
586G4ThreeVector G4Cons::ApproxSurfaceNormal( const G4ThreeVector& p ) const
587{
588 ENorm side ;
589 G4ThreeVector norm ;
590 G4double rho, phi ;
591 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
592 G4double tanRMin, secRMin, pRMin, widRMin ;
593 G4double tanRMax, secRMax, pRMax, widRMax ;
594
595 distZ = std::fabs(std::fabs(p.z()) - fDz) ;
596 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
597
598 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
599 secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
600 pRMin = rho - p.z()*tanRMin ;
601 widRMin = fRmin2 - fDz*tanRMin ;
602 distRMin = std::fabs(pRMin - widRMin)/secRMin ;
603
604 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
605 secRMax = std::sqrt(1+tanRMax*tanRMax) ;
606 pRMax = rho - p.z()*tanRMax ;
607 widRMax = fRmax2 - fDz*tanRMax ;
608 distRMax = std::fabs(pRMax - widRMax)/secRMax ;
609
610 if (distRMin < distRMax) // First minimum
611 {
612 if (distZ < distRMin)
613 {
614 distMin = distZ ;
615 side = kNZ ;
616 }
617 else
618 {
619 distMin = distRMin ;
620 side = kNRMin ;
621 }
622 }
623 else
624 {
625 if (distZ < distRMax)
626 {
627 distMin = distZ ;
628 side = kNZ ;
629 }
630 else
631 {
632 distMin = distRMax ;
633 side = kNRMax ;
634 }
635 }
636 if ( !fPhiFullCone && rho ) // Protected against (0,0,z)
637 {
638 phi = std::atan2(p.y(),p.x()) ;
639
640 if (phi < 0) { phi += twopi; }
641
642 if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
643 else { distSPhi = std::fabs(phi - fSPhi)*rho; }
644
645 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
646
647 // Find new minimum
648
649 if (distSPhi < distEPhi)
650 {
651 if (distSPhi < distMin) { side = kNSPhi; }
652 }
653 else
654 {
655 if (distEPhi < distMin) { side = kNEPhi; }
656 }
657 }
658 switch (side)
659 {
660 case kNRMin: // Inner radius
661 rho *= secRMin ;
662 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
663 break ;
664 case kNRMax: // Outer radius
665 rho *= secRMax ;
666 norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
667 break ;
668 case kNZ: // +/- dz
669 if (p.z() > 0) { norm = G4ThreeVector(0,0,1); }
670 else { norm = G4ThreeVector(0,0,-1); }
671 break ;
672 case kNSPhi:
673 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
674 break ;
675 case kNEPhi:
676 norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
677 break ;
678 default: // Should never reach this case...
679 DumpInfo();
680 G4Exception("G4Cons::ApproxSurfaceNormal()",
681 "GeomSolids1002", JustWarning,
682 "Undefined side for valid surface normal to solid.");
683 break ;
684 }
685 return norm ;
686}
687
688////////////////////////////////////////////////////////////////////////
689//
690// Calculate distance to shape from outside, along normalised vector
691// - return kInfinity if no intersection, or intersection distance <= tolerance
692//
693// - Compute the intersection with the z planes
694// - if at valid r, phi, return
695//
696// -> If point is outside cone, compute intersection with rmax1*0.5
697// - if at valid phi,z return
698// - if inside outer cone, handle case when on tolerant outer cone
699// boundary and heading inwards(->0 to in)
700//
701// -> Compute intersection with inner cone, taking largest +ve root
702// - if valid (in z,phi), save intersction
703//
704// -> If phi segmented, compute intersections with phi half planes
705// - return smallest of valid phi intersections and
706// inner radius intersection
707//
708// NOTE:
709// - `if valid' implies tolerant checking of intersection points
710// - z, phi intersection from Tubs
711
712G4double G4Cons::DistanceToIn( const G4ThreeVector& p,
713 const G4ThreeVector& v ) const
714{
715 G4double snxt = kInfinity ; // snxt = default return value
716 const G4double dRmax = 100*std::min(fRmax1,fRmax2);
717 static const G4double halfCarTolerance=kCarTolerance*0.5;
718 static const G4double halfRadTolerance=kRadTolerance*0.5;
719
720 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; // Data for cones
721 G4double tanRMin,secRMin,rMinAv,rMinOAv ;
722 G4double rout,rin ;
723
724 G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ; // `generous' radii squared
725 G4double tolORMax2,tolIRMax,tolIRMax2 ;
726 G4double tolODz,tolIDz ;
727
728 G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars
729
730 G4double t1,t2,t3,b,c,d ; // Quadratic solver variables
731 G4double nt1,nt2,nt3 ;
732 G4double Comp ;
733
734 G4ThreeVector Normal;
735
736 // Cone Precalcs
737
738 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
739 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
740 rMinAv = (fRmin1 + fRmin2)*0.5 ;
741
742 if (rMinAv > halfRadTolerance)
743 {
744 rMinOAv = rMinAv - halfRadTolerance ;
745 }
746 else
747 {
748 rMinOAv = 0.0 ;
749 }
750 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
751 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
752 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
753 rMaxOAv = rMaxAv + halfRadTolerance ;
754
755 // Intersection with z-surfaces
756
757 tolIDz = fDz - halfCarTolerance ;
758 tolODz = fDz + halfCarTolerance ;
759
760 if (std::fabs(p.z()) >= tolIDz)
761 {
762 if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa
763 {
764 sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance
765
766 if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero
767
768 xi = p.x() + sd*v.x() ; // Intersection coords
769 yi = p.y() + sd*v.y() ;
770 rhoi2 = xi*xi + yi*yi ;
771
772 // Check validity of intersection
773 // Calculate (outer) tolerant radi^2 at intersecion
774
775 if (v.z() > 0)
776 {
777 tolORMin = fRmin1 - halfRadTolerance*secRMin ;
778 tolIRMin = fRmin1 + halfRadTolerance*secRMin ;
779 tolIRMax = fRmax1 - halfRadTolerance*secRMin ;
780 tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
781 (fRmax1 + halfRadTolerance*secRMax) ;
782 }
783 else
784 {
785 tolORMin = fRmin2 - halfRadTolerance*secRMin ;
786 tolIRMin = fRmin2 + halfRadTolerance*secRMin ;
787 tolIRMax = fRmax2 - halfRadTolerance*secRMin ;
788 tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
789 (fRmax2 + halfRadTolerance*secRMax) ;
790 }
791 if ( tolORMin > 0 )
792 {
793 tolORMin2 = tolORMin*tolORMin ;
794 tolIRMin2 = tolIRMin*tolIRMin ;
795 }
796 else
797 {
798 tolORMin2 = 0.0 ;
799 tolIRMin2 = 0.0 ;
800 }
801 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
802 else { tolIRMax2 = 0.0; }
803
804 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
805 {
806 if ( !fPhiFullCone && rhoi2 )
807 {
808 // Psi = angle made with central (average) phi of shape
809
810 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
811
812 if (cosPsi >= cosHDPhiIT) { return sd; }
813 }
814 else
815 {
816 return sd;
817 }
818 }
819 }
820 else // On/outside extent, and heading away -> cannot intersect
821 {
822 return snxt ;
823 }
824 }
825
826// ----> Can not intersect z surfaces
827
828
829// Intersection with outer cone (possible return) and
830// inner cone (must also check phi)
831//
832// Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
833//
834// Intersects with x^2+y^2=(a*z+b)^2
835//
836// where a=tanRMax or tanRMin
837// b=rMaxAv or rMinAv
838//
839// (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
840// t1 t2 t3
841//
842// \--------u-------/ \-----------v----------/ \---------w--------/
843//
844
845 t1 = 1.0 - v.z()*v.z() ;
846 t2 = p.x()*v.x() + p.y()*v.y() ;
847 t3 = p.x()*p.x() + p.y()*p.y() ;
848 rin = tanRMin*p.z() + rMinAv ;
849 rout = tanRMax*p.z() + rMaxAv ;
850
851 // Outer Cone Intersection
852 // Must be outside/on outer cone for valid intersection
853
854 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
855 nt2 = t2 - tanRMax*v.z()*rout ;
856 nt3 = t3 - rout*rout ;
857
858 if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots
859 {
860 b = nt2/nt1;
861 c = nt3/nt1;
862 d = b*b-c ;
863 if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
864 || (rout < 0) )
865 {
866 // If outside real cone (should be rho-rout>kRadTolerance*0.5
867 // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
868
869 if (d >= 0)
870 {
871
872 if ((rout < 0) && (nt3 <= 0))
873 {
874 // Inside `shadow cone' with -ve radius
875 // -> 2nd root could be on real cone
876
877 if (b>0) { sd = c/(-b-std::sqrt(d)); }
878 else { sd = -b + std::sqrt(d); }
879 }
880 else
881 {
882 if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
883 {
884 sd=c/(-b+std::sqrt(d));
885 }
886 else
887 {
888 if ( c <= 0 ) // second >=0
889 {
890 sd = -b + std::sqrt(d) ;
891 }
892 else // both negative, travel away
893 {
894 return kInfinity ;
895 }
896 }
897 }
898 if ( sd > 0 ) // If 'forwards'. Check z intersection
899 {
900 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
901 { // 64 bits systems. Split long distances and recompute
902 G4double fTerm = sd-std::fmod(sd,dRmax);
903 sd = fTerm + DistanceToIn(p+fTerm*v,v);
904 }
905 zi = p.z() + sd*v.z() ;
906
907 if (std::fabs(zi) <= tolODz)
908 {
909 // Z ok. Check phi intersection if reqd
910
911 if ( fPhiFullCone ) { return sd; }
912 else
913 {
914 xi = p.x() + sd*v.x() ;
915 yi = p.y() + sd*v.y() ;
916 ri = rMaxAv + zi*tanRMax ;
917 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
918
919 if ( cosPsi >= cosHDPhiIT ) { return sd; }
920 }
921 }
922 } // end if (sd>0)
923 }
924 }
925 else
926 {
927 // Inside outer cone
928 // check not inside, and heading through G4Cons (-> 0 to in)
929
930 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
931 (rin + halfRadTolerance*secRMin) )
932 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
933 {
934 // Inside cones, delta r -ve, inside z extent
935 // Point is on the Surface => check Direction using Normal.dot(v)
936
937 xi = p.x() ;
938 yi = p.y() ;
939 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
940 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
941 if ( !fPhiFullCone )
942 {
943 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
944 if ( cosPsi >= cosHDPhiIT )
945 {
946 if ( Normal.dot(v) <= 0 ) { return 0.0; }
947 }
948 }
949 else
950 {
951 if ( Normal.dot(v) <= 0 ) { return 0.0; }
952 }
953 }
954 }
955 }
956 else // Single root case
957 {
958 if ( std::fabs(nt2) > kRadTolerance )
959 {
960 sd = -0.5*nt3/nt2 ;
961
962 if ( sd < 0 ) { return kInfinity; } // travel away
963 else // sd >= 0, If 'forwards'. Check z intersection
964 {
965 zi = p.z() + sd*v.z() ;
966
967 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
968 {
969 // Z ok. Check phi intersection if reqd
970
971 if ( fPhiFullCone ) { return sd; }
972 else
973 {
974 xi = p.x() + sd*v.x() ;
975 yi = p.y() + sd*v.y() ;
976 ri = rMaxAv + zi*tanRMax ;
977 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
978
979 if (cosPsi >= cosHDPhiIT) { return sd; }
980 }
981 }
982 }
983 }
984 else // travel || cone surface from its origin
985 {
986 sd = kInfinity ;
987 }
988 }
989
990 // Inner Cone Intersection
991 // o Space is divided into 3 areas:
992 // 1) Radius greater than real inner cone & imaginary cone & outside
993 // tolerance
994 // 2) Radius less than inner or imaginary cone & outside tolarance
995 // 3) Within tolerance of real or imaginary cones
996 // - Extra checks needed for 3's intersections
997 // => lots of duplicated code
998
999 if (rMinAv)
1000 {
1001 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1002 nt2 = t2 - tanRMin*v.z()*rin ;
1003 nt3 = t3 - rin*rin ;
1004
1005 if ( nt1 )
1006 {
1007 if ( nt3 > rin*kRadTolerance*secRMin )
1008 {
1009 // At radius greater than real & imaginary cones
1010 // -> 2nd root, with zi check
1011
1012 b = nt2/nt1 ;
1013 c = nt3/nt1 ;
1014 d = b*b-c ;
1015 if (d >= 0) // > 0
1016 {
1017 if(b>0){sd = c/( -b-std::sqrt(d));}
1018 else {sd = -b + std::sqrt(d) ;}
1019
1020 if ( sd >= 0 ) // > 0
1021 {
1022 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
1023 { // 64 bits systems. Split long distance and recompute
1024 G4double fTerm = sd-std::fmod(sd,dRmax);
1025 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1026 }
1027 zi = p.z() + sd*v.z() ;
1028
1029 if ( std::fabs(zi) <= tolODz )
1030 {
1031 if ( !fPhiFullCone )
1032 {
1033 xi = p.x() + sd*v.x() ;
1034 yi = p.y() + sd*v.y() ;
1035 ri = rMinAv + zi*tanRMin ;
1036 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1037
1038 if (cosPsi >= cosHDPhiIT)
1039 {
1040 if ( sd > halfRadTolerance ) { snxt=sd; }
1041 else
1042 {
1043 // Calculate a normal vector in order to check Direction
1044
1045 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1046 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1047 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1048 }
1049 }
1050 }
1051 else
1052 {
1053 if ( sd > halfRadTolerance ) { return sd; }
1054 else
1055 {
1056 // Calculate a normal vector in order to check Direction
1057
1058 xi = p.x() + sd*v.x() ;
1059 yi = p.y() + sd*v.y() ;
1060 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1061 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1062 if ( Normal.dot(v) <= 0 ) { return sd; }
1063 }
1064 }
1065 }
1066 }
1067 }
1068 }
1069 else if ( nt3 < -rin*kRadTolerance*secRMin )
1070 {
1071 // Within radius of inner cone (real or imaginary)
1072 // -> Try 2nd root, with checking intersection is with real cone
1073 // -> If check fails, try 1st root, also checking intersection is
1074 // on real cone
1075
1076 b = nt2/nt1 ;
1077 c = nt3/nt1 ;
1078 d = b*b - c ;
1079
1080 if ( d >= 0 ) // > 0
1081 {
1082 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1083 else { sd = -b + std::sqrt(d); }
1084 zi = p.z() + sd*v.z() ;
1085 ri = rMinAv + zi*tanRMin ;
1086
1087 if ( ri > 0 )
1088 {
1089 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1090 {
1091 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1092 { // seen on 64 bits systems. Split and recompute
1093 G4double fTerm = sd-std::fmod(sd,dRmax);
1094 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1095 }
1096 if ( !fPhiFullCone )
1097 {
1098 xi = p.x() + sd*v.x() ;
1099 yi = p.y() + sd*v.y() ;
1100 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1101
1102 if (cosPsi >= cosHDPhiOT)
1103 {
1104 if ( sd > halfRadTolerance ) { snxt=sd; }
1105 else
1106 {
1107 // Calculate a normal vector in order to check Direction
1108
1109 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1110 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1111 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1112 }
1113 }
1114 }
1115 else
1116 {
1117 if( sd > halfRadTolerance ) { return sd; }
1118 else
1119 {
1120 // Calculate a normal vector in order to check Direction
1121
1122 xi = p.x() + sd*v.x() ;
1123 yi = p.y() + sd*v.y() ;
1124 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1125 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1126 if ( Normal.dot(v) <= 0 ) { return sd; }
1127 }
1128 }
1129 }
1130 }
1131 else
1132 {
1133 if (b>0) { sd = -b - std::sqrt(d); }
1134 else { sd = c/(-b+std::sqrt(d)); }
1135 zi = p.z() + sd*v.z() ;
1136 ri = rMinAv + zi*tanRMin ;
1137
1138 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1139 {
1140 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1141 { // seen on 64 bits systems. Split and recompute
1142 G4double fTerm = sd-std::fmod(sd,dRmax);
1143 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1144 }
1145 if ( !fPhiFullCone )
1146 {
1147 xi = p.x() + sd*v.x() ;
1148 yi = p.y() + sd*v.y() ;
1149 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1150
1151 if (cosPsi >= cosHDPhiIT)
1152 {
1153 if ( sd > halfRadTolerance ) { snxt=sd; }
1154 else
1155 {
1156 // Calculate a normal vector in order to check Direction
1157
1158 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1159 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1160 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1161 }
1162 }
1163 }
1164 else
1165 {
1166 if ( sd > halfRadTolerance ) { return sd; }
1167 else
1168 {
1169 // Calculate a normal vector in order to check Direction
1170
1171 xi = p.x() + sd*v.x() ;
1172 yi = p.y() + sd*v.y() ;
1173 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1174 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1175 if ( Normal.dot(v) <= 0 ) { return sd; }
1176 }
1177 }
1178 }
1179 }
1180 }
1181 }
1182 else
1183 {
1184 // Within kRadTol*0.5 of inner cone (real OR imaginary)
1185 // ----> Check not travelling through (=>0 to in)
1186 // ----> if not:
1187 // -2nd root with validity check
1188
1189 if ( std::fabs(p.z()) <= tolODz )
1190 {
1191 if ( nt2 > 0 )
1192 {
1193 // Inside inner real cone, heading outwards, inside z range
1194
1195 if ( !fPhiFullCone )
1196 {
1197 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
1198
1199 if (cosPsi >= cosHDPhiIT) { return 0.0; }
1200 }
1201 else { return 0.0; }
1202 }
1203 else
1204 {
1205 // Within z extent, but not travelling through
1206 // -> 2nd root or kInfinity if 1st root on imaginary cone
1207
1208 b = nt2/nt1 ;
1209 c = nt3/nt1 ;
1210 d = b*b - c ;
1211
1212 if ( d >= 0 ) // > 0
1213 {
1214 if (b>0) { sd = -b - std::sqrt(d); }
1215 else { sd = c/(-b+std::sqrt(d)); }
1216 zi = p.z() + sd*v.z() ;
1217 ri = rMinAv + zi*tanRMin ;
1218
1219 if ( ri > 0 ) // 2nd root
1220 {
1221 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1222 else { sd = -b + std::sqrt(d); }
1223
1224 zi = p.z() + sd*v.z() ;
1225
1226 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1227 {
1228 if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1229 { // seen on 64 bits systems. Split and recompute
1230 G4double fTerm = sd-std::fmod(sd,dRmax);
1231 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1232 }
1233 if ( !fPhiFullCone )
1234 {
1235 xi = p.x() + sd*v.x() ;
1236 yi = p.y() + sd*v.y() ;
1237 ri = rMinAv + zi*tanRMin ;
1238 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1239
1240 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1241 }
1242 else { return sd; }
1243 }
1244 }
1245 else { return kInfinity; }
1246 }
1247 }
1248 }
1249 else // 2nd root
1250 {
1251 b = nt2/nt1 ;
1252 c = nt3/nt1 ;
1253 d = b*b - c ;
1254
1255 if ( d > 0 )
1256 {
1257 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1258 else { sd = -b + std::sqrt(d) ; }
1259 zi = p.z() + sd*v.z() ;
1260
1261 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1262 {
1263 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1264 { // seen on 64 bits systems. Split and recompute
1265 G4double fTerm = sd-std::fmod(sd,dRmax);
1266 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1267 }
1268 if ( !fPhiFullCone )
1269 {
1270 xi = p.x() + sd*v.x();
1271 yi = p.y() + sd*v.y();
1272 ri = rMinAv + zi*tanRMin ;
1273 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1274
1275 if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1276 }
1277 else { return sd; }
1278 }
1279 }
1280 }
1281 }
1282 }
1283 }
1284
1285 // Phi segment intersection
1286 //
1287 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1288 //
1289 // o NOTE: Large duplication of code between sphi & ephi checks
1290 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1291 // intersection check <=0 -> >=0
1292 // -> Should use some form of loop Construct
1293
1294 if ( !fPhiFullCone )
1295 {
1296 // First phi surface (starting phi)
1297
1298 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1299
1300 if ( Comp < 0 ) // Component in outwards normal dirn
1301 {
1302 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1303
1304 if (Dist < halfCarTolerance)
1305 {
1306 sd = Dist/Comp ;
1307
1308 if ( sd < snxt )
1309 {
1310 if ( sd < 0 ) { sd = 0.0; }
1311
1312 zi = p.z() + sd*v.z() ;
1313
1314 if ( std::fabs(zi) <= tolODz )
1315 {
1316 xi = p.x() + sd*v.x() ;
1317 yi = p.y() + sd*v.y() ;
1318 rhoi2 = xi*xi + yi*yi ;
1319 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1320 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1321
1322 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1323 {
1324 // z and r intersections good - check intersecting with
1325 // correct half-plane
1326
1327 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1328 }
1329 }
1330 }
1331 }
1332 }
1333
1334 // Second phi surface (Ending phi)
1335
1336 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1337
1338 if ( Comp < 0 ) // Component in outwards normal dirn
1339 {
1340 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1341 if (Dist < halfCarTolerance)
1342 {
1343 sd = Dist/Comp ;
1344
1345 if ( sd < snxt )
1346 {
1347 if ( sd < 0 ) { sd = 0.0; }
1348
1349 zi = p.z() + sd*v.z() ;
1350
1351 if (std::fabs(zi) <= tolODz)
1352 {
1353 xi = p.x() + sd*v.x() ;
1354 yi = p.y() + sd*v.y() ;
1355 rhoi2 = xi*xi + yi*yi ;
1356 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1357 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1358
1359 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1360 {
1361 // z and r intersections good - check intersecting with
1362 // correct half-plane
1363
1364 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1365 }
1366 }
1367 }
1368 }
1369 }
1370 }
1371 if (snxt < halfCarTolerance) { snxt = 0.; }
1372
1373 return snxt ;
1374}
1375
1376//////////////////////////////////////////////////////////////////////////////
1377//
1378// Calculate distance (<= actual) to closest surface of shape from outside
1379// - Calculate distance to z, radial planes
1380// - Only to phi planes if outside phi extent
1381// - Return 0 if point inside
1382
1383G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
1384{
1385 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1386 G4double tanRMin, secRMin, pRMin ;
1387 G4double tanRMax, secRMax, pRMax ;
1388
1389 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1390 safeZ = std::fabs(p.z()) - fDz ;
1391
1392 if ( fRmin1 || fRmin2 )
1393 {
1394 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1395 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1396 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1397 safeR1 = (pRMin - rho)/secRMin ;
1398
1399 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1400 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1401 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1402 safeR2 = (rho - pRMax)/secRMax ;
1403
1404 if ( safeR1 > safeR2) { safe = safeR1; }
1405 else { safe = safeR2; }
1406 }
1407 else
1408 {
1409 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1410 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1411 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1412 safe = (rho - pRMax)/secRMax ;
1413 }
1414 if ( safeZ > safe ) { safe = safeZ; }
1415
1416 if ( !fPhiFullCone && rho )
1417 {
1418 // Psi=angle from central phi to point
1419
1420 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1421
1422 if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1423 {
1424 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1425 {
1426 safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1427 }
1428 else
1429 {
1430 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1431 }
1432 if ( safePhi > safe ) { safe = safePhi; }
1433 }
1434 }
1435 if ( safe < 0.0 ) { safe = 0.0; }
1436
1437 return safe ;
1438}
1439
1440///////////////////////////////////////////////////////////////
1441//
1442// Calculate distance to surface of shape from 'inside', allowing for tolerance
1443// - Only Calc rmax intersection if no valid rmin intersection
1444
1445G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
1446 const G4ThreeVector& v,
1447 const G4bool calcNorm,
1448 G4bool *validNorm,
1449 G4ThreeVector *n) const
1450{
1451 ESide side = kNull, sider = kNull, sidephi = kNull;
1452
1453 static const G4double halfCarTolerance=kCarTolerance*0.5;
1454 static const G4double halfRadTolerance=kRadTolerance*0.5;
1455 static const G4double halfAngTolerance=kAngTolerance*0.5;
1456
1457 G4double snxt,srd,sphi,pdist ;
1458
1459 G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone
1460 G4double tanRMin, secRMin, rMinAv ; // Data for inner cone
1461
1462 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1463 G4double b, c, d, sr2, sr3 ;
1464
1465 // Vars for intersection within tolerance
1466
1467 ESide sidetol = kNull ;
1468 G4double slentol = kInfinity ;
1469
1470 // Vars for phi intersection:
1471
1472 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1473 G4double zi, ri, deltaRoi2 ;
1474
1475 // Z plane intersection
1476
1477 if ( v.z() > 0.0 )
1478 {
1479 pdist = fDz - p.z() ;
1480
1481 if (pdist > halfCarTolerance)
1482 {
1483 snxt = pdist/v.z() ;
1484 side = kPZ ;
1485 }
1486 else
1487 {
1488 if (calcNorm)
1489 {
1490 *n = G4ThreeVector(0,0,1) ;
1491 *validNorm = true ;
1492 }
1493 return snxt = 0.0;
1494 }
1495 }
1496 else if ( v.z() < 0.0 )
1497 {
1498 pdist = fDz + p.z() ;
1499
1500 if ( pdist > halfCarTolerance)
1501 {
1502 snxt = -pdist/v.z() ;
1503 side = kMZ ;
1504 }
1505 else
1506 {
1507 if ( calcNorm )
1508 {
1509 *n = G4ThreeVector(0,0,-1) ;
1510 *validNorm = true ;
1511 }
1512 return snxt = 0.0 ;
1513 }
1514 }
1515 else // Travel perpendicular to z axis
1516 {
1517 snxt = kInfinity ;
1518 side = kNull ;
1519 }
1520
1521 // Radial Intersections
1522 //
1523 // Intersection with outer cone (possible return) and
1524 // inner cone (must also check phi)
1525 //
1526 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1527 //
1528 // Intersects with x^2+y^2=(a*z+b)^2
1529 //
1530 // where a=tanRMax or tanRMin
1531 // b=rMaxAv or rMinAv
1532 //
1533 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
1534 // t1 t2 t3
1535 //
1536 // \--------u-------/ \-----------v----------/ \---------w--------/
1537
1538 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1539 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1540 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
1541
1542
1543 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1544 t2 = p.x()*v.x() + p.y()*v.y() ;
1545 t3 = p.x()*p.x() + p.y()*p.y() ;
1546 rout = tanRMax*p.z() + rMaxAv ;
1547
1548 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1549 nt2 = t2 - tanRMax*v.z()*rout ;
1550 nt3 = t3 - rout*rout ;
1551
1552 if (v.z() > 0.0)
1553 {
1554 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1555 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1556 }
1557 else if ( v.z() < 0.0 )
1558 {
1559 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1560 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1561 }
1562 else
1563 {
1564 deltaRoi2 = 1.0;
1565 }
1566
1567 if ( nt1 && (deltaRoi2 > 0.0) )
1568 {
1569 // Equation quadratic => 2 roots : second root must be leaving
1570
1571 b = nt2/nt1 ;
1572 c = nt3/nt1 ;
1573 d = b*b - c ;
1574
1575 if ( d >= 0 )
1576 {
1577 // Check if on outer cone & heading outwards
1578 // NOTE: Should use rho-rout>-kRadTolerance*0.5
1579
1580 if (nt3 > -halfRadTolerance && nt2 >= 0 )
1581 {
1582 if (calcNorm)
1583 {
1584 risec = std::sqrt(t3)*secRMax ;
1585 *validNorm = true ;
1586 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1587 }
1588 return snxt=0 ;
1589 }
1590 else
1591 {
1592 sider = kRMax ;
1593 if (b>0) { srd = -b - std::sqrt(d); }
1594 else { srd = c/(-b+std::sqrt(d)) ; }
1595
1596 zi = p.z() + srd*v.z() ;
1597 ri = tanRMax*zi + rMaxAv ;
1598
1599 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1600 {
1601 // An intersection within the tolerance
1602 // we will Store it in case it is good -
1603 //
1604 slentol = srd ;
1605 sidetol = kRMax ;
1606 }
1607 if ( (ri < 0) || (srd < halfRadTolerance) )
1608 {
1609 // Safety: if both roots -ve ensure that srd cannot `win'
1610 // distance to out
1611
1612 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1613 else { sr2 = -b + std::sqrt(d); }
1614 zi = p.z() + sr2*v.z() ;
1615 ri = tanRMax*zi + rMaxAv ;
1616
1617 if ((ri >= 0) && (sr2 > halfRadTolerance))
1618 {
1619 srd = sr2;
1620 }
1621 else
1622 {
1623 srd = kInfinity ;
1624
1625 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
1626 {
1627 // An intersection within the tolerance.
1628 // Storing it in case it is good.
1629
1630 slentol = sr2 ;
1631 sidetol = kRMax ;
1632 }
1633 }
1634 }
1635 }
1636 }
1637 else
1638 {
1639 // No intersection with outer cone & not parallel
1640 // -> already outside, no intersection
1641
1642 if ( calcNorm )
1643 {
1644 risec = std::sqrt(t3)*secRMax;
1645 *validNorm = true;
1646 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1647 }
1648 return snxt = 0.0 ;
1649 }
1650 }
1651 else if ( nt2 && (deltaRoi2 > 0.0) )
1652 {
1653 // Linear case (only one intersection) => point outside outer cone
1654
1655 if ( calcNorm )
1656 {
1657 risec = std::sqrt(t3)*secRMax;
1658 *validNorm = true;
1659 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1660 }
1661 return snxt = 0.0 ;
1662 }
1663 else
1664 {
1665 // No intersection -> parallel to outer cone
1666 // => Z or inner cone intersection
1667
1668 srd = kInfinity ;
1669 }
1670
1671 // Check possible intersection within tolerance
1672
1673 if ( slentol <= halfCarTolerance )
1674 {
1675 // An intersection within the tolerance was found.
1676 // We must accept it only if the momentum points outwards.
1677 //
1678 // G4ThreeVector ptTol ; // The point of the intersection
1679 // ptTol= p + slentol*v ;
1680 // ri=tanRMax*zi+rMaxAv ;
1681 //
1682 // Calculate a normal vector, as below
1683
1684 xi = p.x() + slentol*v.x();
1685 yi = p.y() + slentol*v.y();
1686 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1687 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
1688
1689 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly
1690 {
1691 if ( calcNorm )
1692 {
1693 *n = Normal.unit() ;
1694 *validNorm = true ;
1695 }
1696 return snxt = 0.0 ;
1697 }
1698 else // On the surface, but not heading out so we ignore this intersection
1699 { // (as it is within tolerance).
1700 slentol = kInfinity ;
1701 }
1702 }
1703
1704 // Inner Cone intersection
1705
1706 if ( fRmin1 || fRmin2 )
1707 {
1708 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1709 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1710
1711 if ( nt1 )
1712 {
1713 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1714 rMinAv = (fRmin1 + fRmin2)*0.5 ;
1715 rin = tanRMin*p.z() + rMinAv ;
1716 nt2 = t2 - tanRMin*v.z()*rin ;
1717 nt3 = t3 - rin*rin ;
1718
1719 // Equation quadratic => 2 roots : first root must be leaving
1720
1721 b = nt2/nt1 ;
1722 c = nt3/nt1 ;
1723 d = b*b - c ;
1724
1725 if ( d >= 0.0 )
1726 {
1727 // NOTE: should be rho-rin<kRadTolerance*0.5,
1728 // but using squared versions for efficiency
1729
1730 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
1731 {
1732 if ( nt2 < 0.0 )
1733 {
1734 if (calcNorm) { *validNorm = false; }
1735 return snxt = 0.0;
1736 }
1737 }
1738 else
1739 {
1740 if (b>0) { sr2 = -b - std::sqrt(d); }
1741 else { sr2 = c/(-b+std::sqrt(d)); }
1742 zi = p.z() + sr2*v.z() ;
1743 ri = tanRMin*zi + rMinAv ;
1744
1745 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
1746 {
1747 // An intersection within the tolerance
1748 // storing it in case it is good.
1749
1750 slentol = sr2 ;
1751 sidetol = kRMax ;
1752 }
1753 if( (ri<0) || (sr2 < halfRadTolerance) )
1754 {
1755 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1756 else { sr3 = -b + std::sqrt(d) ; }
1757
1758 // Safety: if both roots -ve ensure that srd cannot `win'
1759 // distancetoout
1760
1761 if ( sr3 > halfRadTolerance )
1762 {
1763 if( sr3 < srd )
1764 {
1765 zi = p.z() + sr3*v.z() ;
1766 ri = tanRMin*zi + rMinAv ;
1767
1768 if ( ri >= 0.0 )
1769 {
1770 srd=sr3 ;
1771 sider=kRMin ;
1772 }
1773 }
1774 }
1775 else if ( sr3 > -halfRadTolerance )
1776 {
1777 // Intersection in tolerance. Store to check if it's good
1778
1779 slentol = sr3 ;
1780 sidetol = kRMin ;
1781 }
1782 }
1783 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
1784 {
1785 srd = sr2 ;
1786 sider = kRMin ;
1787 }
1788 else if (sr2 > -halfCarTolerance)
1789 {
1790 // Intersection in tolerance. Store to check if it's good
1791
1792 slentol = sr2 ;
1793 sidetol = kRMin ;
1794 }
1795 if( slentol <= halfCarTolerance )
1796 {
1797 // An intersection within the tolerance was found.
1798 // We must accept it only if the momentum points outwards.
1799
1800 G4ThreeVector Normal ;
1801
1802 // Calculate a normal vector, as below
1803
1804 xi = p.x() + slentol*v.x() ;
1805 yi = p.y() + slentol*v.y() ;
1806 if( sidetol==kRMax )
1807 {
1808 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1809 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1810 }
1811 else
1812 {
1813 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1814 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1815 }
1816 if( Normal.dot(v) > 0 )
1817 {
1818 // We will leave the cone immediately
1819
1820 if( calcNorm )
1821 {
1822 *n = Normal.unit() ;
1823 *validNorm = true ;
1824 }
1825 return snxt = 0.0 ;
1826 }
1827 else
1828 {
1829 // On the surface, but not heading out so we ignore this
1830 // intersection (as it is within tolerance).
1831
1832 slentol = kInfinity ;
1833 }
1834 }
1835 }
1836 }
1837 }
1838 }
1839
1840 // Linear case => point outside inner cone ---> outer cone intersect
1841 //
1842 // Phi Intersection
1843
1844 if ( !fPhiFullCone )
1845 {
1846 // add angle calculation with correction
1847 // of the difference in domain of atan2 and Sphi
1848
1849 vphi = std::atan2(v.y(),v.x()) ;
1850
1851 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1852 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1853
1854 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1855 {
1856 // pDist -ve when inside
1857
1858 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1859 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1860
1861 // Comp -ve when in direction of outwards normal
1862
1863 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1864 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1865
1866 sidephi = kNull ;
1867
1868 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1869 && (pDistE <= halfCarTolerance) ) )
1870 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1871 && (pDistE > halfCarTolerance) ) ) )
1872 {
1873 // Inside both phi *full* planes
1874 if ( compS < 0 )
1875 {
1876 sphi = pDistS/compS ;
1877 if (sphi >= -halfCarTolerance)
1878 {
1879 xi = p.x() + sphi*v.x() ;
1880 yi = p.y() + sphi*v.y() ;
1881
1882 // Check intersecting with correct half-plane
1883 // (if not -> no intersect)
1884 //
1885 if ( (std::fabs(xi)<=kCarTolerance)
1886 && (std::fabs(yi)<=kCarTolerance) )
1887 {
1888 sidephi= kSPhi;
1889 if ( ( fSPhi-halfAngTolerance <= vphi )
1890 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
1891 {
1892 sphi = kInfinity;
1893 }
1894 }
1895 else
1896 if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
1897 {
1898 sphi = kInfinity ;
1899 }
1900 else
1901 {
1902 sidephi = kSPhi ;
1903 if ( pDistS > -halfCarTolerance )
1904 {
1905 sphi = 0.0 ; // Leave by sphi immediately
1906 }
1907 }
1908 }
1909 else
1910 {
1911 sphi = kInfinity ;
1912 }
1913 }
1914 else
1915 {
1916 sphi = kInfinity ;
1917 }
1918
1919 if ( compE < 0 )
1920 {
1921 sphi2 = pDistE/compE ;
1922
1923 // Only check further if < starting phi intersection
1924 //
1925 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1926 {
1927 xi = p.x() + sphi2*v.x() ;
1928 yi = p.y() + sphi2*v.y() ;
1929
1930 // Check intersecting with correct half-plane
1931
1932 if ( (std::fabs(xi)<=kCarTolerance)
1933 && (std::fabs(yi)<=kCarTolerance) )
1934 {
1935 // Leaving via ending phi
1936
1937 if(!( (fSPhi-halfAngTolerance <= vphi)
1938 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
1939 {
1940 sidephi = kEPhi ;
1941 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1942 else { sphi = 0.0; }
1943 }
1944 }
1945 else // Check intersecting with correct half-plane
1946 if ( yi*cosCPhi-xi*sinCPhi >= 0 )
1947 {
1948 // Leaving via ending phi
1949
1950 sidephi = kEPhi ;
1951 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1952 else { sphi = 0.0; }
1953 }
1954 }
1955 }
1956 }
1957 else
1958 {
1959 sphi = kInfinity ;
1960 }
1961 }
1962 else
1963 {
1964 // On z axis + travel not || to z axis -> if phi of vector direction
1965 // within phi of shape, Step limited by rmax, else Step =0
1966
1967 if ( (fSPhi-halfAngTolerance <= vphi)
1968 && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
1969 {
1970 sphi = kInfinity ;
1971 }
1972 else
1973 {
1974 sidephi = kSPhi ; // arbitrary
1975 sphi = 0.0 ;
1976 }
1977 }
1978 if ( sphi < snxt ) // Order intersecttions
1979 {
1980 snxt=sphi ;
1981 side=sidephi ;
1982 }
1983 }
1984 if ( srd < snxt ) // Order intersections
1985 {
1986 snxt = srd ;
1987 side = sider ;
1988 }
1989 if (calcNorm)
1990 {
1991 switch(side)
1992 { // Note: returned vector not normalised
1993 case kRMax: // (divide by frmax for unit vector)
1994 xi = p.x() + snxt*v.x() ;
1995 yi = p.y() + snxt*v.y() ;
1996 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1997 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1998 *validNorm = true ;
1999 break ;
2000 case kRMin:
2001 *validNorm = false ; // Rmin is inconvex
2002 break ;
2003 case kSPhi:
2004 if ( fDPhi <= pi )
2005 {
2006 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0);
2007 *validNorm = true ;
2008 }
2009 else
2010 {
2011 *validNorm = false ;
2012 }
2013 break ;
2014 case kEPhi:
2015 if ( fDPhi <= pi )
2016 {
2017 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
2018 *validNorm = true ;
2019 }
2020 else
2021 {
2022 *validNorm = false ;
2023 }
2024 break ;
2025 case kPZ:
2026 *n = G4ThreeVector(0,0,1) ;
2027 *validNorm = true ;
2028 break ;
2029 case kMZ:
2030 *n = G4ThreeVector(0,0,-1) ;
2031 *validNorm = true ;
2032 break ;
2033 default:
2034 G4cout << G4endl ;
2035 DumpInfo();
2036 std::ostringstream message;
2037 G4int oldprc = message.precision(16) ;
2038 message << "Undefined side for valid surface normal to solid."
2039 << G4endl
2040 << "Position:" << G4endl << G4endl
2041 << "p.x() = " << p.x()/mm << " mm" << G4endl
2042 << "p.y() = " << p.y()/mm << " mm" << G4endl
2043 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2044 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2045 << " mm" << G4endl << G4endl ;
2046 if( p.x() != 0. || p.x() != 0.)
2047 {
2048 message << "point phi = " << std::atan2(p.y(),p.x())/degree
2049 << " degree" << G4endl << G4endl ;
2050 }
2051 message << "Direction:" << G4endl << G4endl
2052 << "v.x() = " << v.x() << G4endl
2053 << "v.y() = " << v.y() << G4endl
2054 << "v.z() = " << v.z() << G4endl<< G4endl
2055 << "Proposed distance :" << G4endl<< G4endl
2056 << "snxt = " << snxt/mm << " mm" << G4endl ;
2057 message.precision(oldprc) ;
2058 G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002",
2059 JustWarning, message) ;
2060 break ;
2061 }
2062 }
2063 if (snxt < halfCarTolerance) { snxt = 0.; }
2064
2065 return snxt ;
2066}
2067
2068//////////////////////////////////////////////////////////////////
2069//
2070// Calculate distance (<=actual) to closest surface of shape from inside
2071
2072G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const
2073{
2074 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2075 G4double tanRMin, secRMin, pRMin;
2076 G4double tanRMax, secRMax, pRMax;
2077
2078#ifdef G4CSGDEBUG
2079 if( Inside(p) == kOutside )
2080 {
2081 G4int oldprc=G4cout.precision(16) ;
2082 G4cout << G4endl ;
2083 DumpInfo();
2084 G4cout << "Position:" << G4endl << G4endl ;
2085 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2086 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2087 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2088 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2089 << " mm" << G4endl << G4endl ;
2090 if( (p.x() != 0.) || (p.x() != 0.) )
2091 {
2092 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2093 << " degree" << G4endl << G4endl ;
2094 }
2095 G4cout.precision(oldprc) ;
2096 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
2097 JustWarning, "Point p is outside !?" );
2098 }
2099#endif
2100
2101 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2102 safeZ = fDz - std::fabs(p.z()) ;
2103
2104 if (fRmin1 || fRmin2)
2105 {
2106 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2107 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2108 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2109 safeR1 = (rho - pRMin)/secRMin ;
2110 }
2111 else
2112 {
2113 safeR1 = kInfinity ;
2114 }
2115
2116 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2117 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2118 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2119 safeR2 = (pRMax - rho)/secRMax ;
2120
2121 if (safeR1 < safeR2) { safe = safeR1; }
2122 else { safe = safeR2; }
2123 if (safeZ < safe) { safe = safeZ ; }
2124
2125 // Check if phi divided, Calc distances closest phi plane
2126
2127 if (!fPhiFullCone)
2128 {
2129 // Above/below central phi of G4Cons?
2130
2131 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2132 {
2133 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2134 }
2135 else
2136 {
2137 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2138 }
2139 if (safePhi < safe) { safe = safePhi; }
2140 }
2141 if ( safe < 0 ) { safe = 0; }
2142
2143 return safe ;
2144}
2145
2146////////////////////////////////////////////////////////////////////////////
2147//
2148// Create a List containing the transformed vertices
2149// Ordering [0-3] -fDz cross section
2150// [4-7] +fDz cross section such that [0] is below [4],
2151// [1] below [5] etc.
2152// Note:
2153// Caller has deletion resposibility
2154// Potential improvement: For last slice, use actual ending angle
2155// to avoid rounding error problems.
2156
2157G4ThreeVectorList*
2158G4Cons::CreateRotatedVertices(const G4AffineTransform& pTransform) const
2159{
2160 G4ThreeVectorList* vertices ;
2161 G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
2162 G4double meshAngle, meshRMax1, meshRMax2, crossAngle;
2163 G4double cosCrossAngle, sinCrossAngle, sAngle ;
2164 G4double rMaxX1, rMaxX2, rMaxY1, rMaxY2, rMinX1, rMinX2, rMinY1, rMinY2 ;
2165 G4int crossSection, noCrossSections ;
2166
2167 // Compute no of cross-sections necessary to mesh cone
2168
2169 noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
2170
2171 if (noCrossSections < kMinMeshSections)
2172 {
2173 noCrossSections = kMinMeshSections ;
2174 }
2175 else if (noCrossSections > kMaxMeshSections)
2176 {
2177 noCrossSections = kMaxMeshSections ;
2178 }
2179 meshAngle = fDPhi/(noCrossSections - 1) ;
2180
2181 meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ;
2182 meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ;
2183
2184 // If complete in phi, set start angle such that mesh will be at RMax
2185 // on the x axis. Will give better extent calculations when not rotated.
2186
2187 if ( fPhiFullCone && (fSPhi == 0.0) )
2188 {
2189 sAngle = -meshAngle*0.5 ;
2190 }
2191 else
2192 {
2193 sAngle = fSPhi ;
2194 }
2195 vertices = new G4ThreeVectorList();
2196
2197 if (vertices)
2198 {
2199 vertices->reserve(noCrossSections*4) ;
2200 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++)
2201 {
2202 // Compute coordinates of cross section at section crossSection
2203
2204 crossAngle = sAngle + crossSection*meshAngle ;
2205 cosCrossAngle = std::cos(crossAngle) ;
2206 sinCrossAngle = std::sin(crossAngle) ;
2207
2208 rMaxX1 = meshRMax1*cosCrossAngle ;
2209 rMaxY1 = meshRMax1*sinCrossAngle ;
2210 rMaxX2 = meshRMax2*cosCrossAngle ;
2211 rMaxY2 = meshRMax2*sinCrossAngle ;
2212
2213 rMinX1 = fRmin1*cosCrossAngle ;
2214 rMinY1 = fRmin1*sinCrossAngle ;
2215 rMinX2 = fRmin2*cosCrossAngle ;
2216 rMinY2 = fRmin2*sinCrossAngle ;
2217
2218 vertex0 = G4ThreeVector(rMinX1,rMinY1,-fDz) ;
2219 vertex1 = G4ThreeVector(rMaxX1,rMaxY1,-fDz) ;
2220 vertex2 = G4ThreeVector(rMaxX2,rMaxY2,+fDz) ;
2221 vertex3 = G4ThreeVector(rMinX2,rMinY2,+fDz) ;
2222
2223 vertices->push_back(pTransform.TransformPoint(vertex0)) ;
2224 vertices->push_back(pTransform.TransformPoint(vertex1)) ;
2225 vertices->push_back(pTransform.TransformPoint(vertex2)) ;
2226 vertices->push_back(pTransform.TransformPoint(vertex3)) ;
2227 }
2228 }
2229 else
2230 {
2231 DumpInfo();
2232 G4Exception("G4Cons::CreateRotatedVertices()",
2233 "GeomSolids0003", FatalException,
2234 "Error in allocation of vertices. Out of memory !");
2235 }
2236
2237 return vertices ;
2238}
2239
2240//////////////////////////////////////////////////////////////////////////
2241//
2242// GetEntityType
2243
2244G4GeometryType G4Cons::GetEntityType() const
2245{
2246 return G4String("G4Cons");
2247}
2248
2249//////////////////////////////////////////////////////////////////////////
2250//
2251// Make a clone of the object
2252//
2253G4VSolid* G4Cons::Clone() const
2254{
2255 return new G4Cons(*this);
2256}
2257
2258//////////////////////////////////////////////////////////////////////////
2259//
2260// Stream object contents to an output stream
2261
2262std::ostream& G4Cons::StreamInfo(std::ostream& os) const
2263{
2264 G4int oldprc = os.precision(16);
2265 os << "-----------------------------------------------------------\n"
2266 << " *** Dump for solid - " << GetName() << " ***\n"
2267 << " ===================================================\n"
2268 << " Solid type: G4Cons\n"
2269 << " Parameters: \n"
2270 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2271 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2272 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2273 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2274 << " half length in Z : " << fDz/mm << " mm \n"
2275 << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2276 << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2277 << "-----------------------------------------------------------\n";
2278 os.precision(oldprc);
2279
2280 return os;
2281}
2282
2283
2284
2285/////////////////////////////////////////////////////////////////////////
2286//
2287// GetPointOnSurface
2288
2289G4ThreeVector G4Cons::GetPointOnSurface() const
2290{
2291 // declare working variables
2292 //
2293 G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2294 G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2295 rone = (fRmax1-fRmax2)/(2.*fDz);
2296 rtwo = (fRmin1-fRmin2)/(2.*fDz);
2297 qone=0.; qtwo=0.;
2298 if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); }
2299 if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); }
2300 slin = std::sqrt(sqr(fRmin1-fRmin2)+sqr(2.*fDz));
2301 slout = std::sqrt(sqr(fRmax1-fRmax2)+sqr(2.*fDz));
2302 Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout;
2303 Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin;
2304 Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1);
2305 Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2);
2306 Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
2307
2308 phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
2309 cosu = std::cos(phi); sinu = std::sin(phi);
2310 rRand1 = GetRadiusInRing(fRmin1, fRmin2);
2311 rRand2 = GetRadiusInRing(fRmax1, fRmax2);
2312
2313 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
2314 chose = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2315
2316 if( (chose >= 0.) && (chose < Aone) )
2317 {
2318 if(fRmin1 != fRmin2)
2319 {
2320 zRand = RandFlat::shoot(-1.*fDz,fDz);
2321 return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
2322 rtwo*sinu*(qtwo-zRand), zRand);
2323 }
2324 else
2325 {
2326 return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
2327 RandFlat::shoot(-1.*fDz,fDz));
2328 }
2329 }
2330 else if( (chose >= Aone) && (chose <= Aone + Atwo) )
2331 {
2332 if(fRmax1 != fRmax2)
2333 {
2334 zRand = RandFlat::shoot(-1.*fDz,fDz);
2335 return G4ThreeVector (rone*cosu*(qone-zRand),
2336 rone*sinu*(qone-zRand), zRand);
2337 }
2338 else
2339 {
2340 return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
2341 RandFlat::shoot(-1.*fDz,fDz));
2342 }
2343 }
2344 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
2345 {
2346 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
2347 }
2348 else if( (chose >= Aone + Atwo + Athree)
2349 && (chose < Aone + Atwo + Athree + Afour) )
2350 {
2351 return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz);
2352 }
2353 else if( (chose >= Aone + Atwo + Athree + Afour)
2354 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2355 {
2356 zRand = RandFlat::shoot(-1.*fDz,fDz);
2357 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2358 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2359 return G4ThreeVector (rRand1*std::cos(fSPhi),
2360 rRand1*std::sin(fSPhi), zRand);
2361 }
2362 else
2363 {
2364 zRand = RandFlat::shoot(-1.*fDz,fDz);
2365 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2366 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2367 return G4ThreeVector (rRand1*std::cos(fSPhi+fDPhi),
2368 rRand1*std::sin(fSPhi+fDPhi), zRand);
2369 }
2370}
2371
2372//////////////////////////////////////////////////////////////////////////
2373//
2374// Methods for visualisation
2375
2376void G4Cons::DescribeYourselfTo (G4VGraphicsScene& scene) const
2377{
2378 scene.AddSolid (*this);
2379}
2380
2381G4Polyhedron* G4Cons::CreatePolyhedron () const
2382{
2383 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
2384}
2385
2386G4NURBS* G4Cons::CreateNURBS () const
2387{
2388 G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ;
2389 return new G4NURBSbox (RMax, RMax, fDz); // Box for now!!!
2390}
ESide
ESide
Definition: G4Cons.cc:68
kPZ
@ kPZ
Definition: G4Cons.cc:68
kRMax
@ kRMax
Definition: G4Cons.cc:68
kNull
@ kNull
Definition: G4Cons.cc:68
kSPhi
@ kSPhi
Definition: G4Cons.cc:68
kRMin
@ kRMin
Definition: G4Cons.cc:68
kEPhi
@ kEPhi
Definition: G4Cons.cc:68
kMZ
@ kMZ
Definition: G4Cons.cc:68
ENorm
ENorm
Definition: G4Cons.cc:72
kNRMax
@ kNRMax
Definition: G4Cons.cc:72
kNEPhi
@ kNEPhi
Definition: G4Cons.cc:72
kNSPhi
@ kNSPhi
Definition: G4Cons.cc:72
kNRMin
@ kNRMin
Definition: G4Cons.cc:72
kNZ
@ kNZ
Definition: G4Cons.cc:72
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::RandFlat::shoot
static double shoot()
Definition: RandFlat.cc:59
G4AffineTransform::IsRotated
G4bool IsRotated() const
G4AffineTransform::Inverse
G4AffineTransform Inverse() const
G4AffineTransform::NetTranslation
G4ThreeVector NetTranslation() const
G4AffineTransform::TransformPoint
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4CSGSolid::GetRadiusInRing
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101
G4CSGSolid::operator=
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
G4Cons
Definition: G4Cons.hh:75
G4Cons::Clone
G4VSolid * Clone() const
Definition: G4Cons.cc:2253
G4Cons::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Cons.cc:249
G4Cons::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Cons.cc:2376
G4Cons::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4Cons.cc:2244
G4Cons::~G4Cons
~G4Cons()
Definition: G4Cons.cc:138
G4Cons::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4Cons.cc:2381
G4Cons::G4Cons
G4Cons(const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4Cons.cc:79
G4Cons::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Cons.cc:479
G4Cons::operator=
G4Cons & operator=(const G4Cons &rhs)
Definition: G4Cons.cc:161
G4Cons::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Cons.cc:1445
G4Cons::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4Cons.cc:191
G4Cons::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4Cons.cc:2289
G4Cons::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Cons.cc:261
G4Cons::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Cons.cc:712
G4Cons::CreateNURBS
G4NURBS * CreateNURBS() const
Definition: G4Cons.cc:2386
G4Cons::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Cons.cc:2262
G4GeometryTolerance::GetRadialTolerance
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
G4GeometryTolerance::GetAngularTolerance
G4double GetAngularTolerance() const
Definition: G4GeometryTolerance.cc:92
G4VGraphicsScene::AddSolid
virtual void AddSolid(const G4Box &)=0
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:94
G4VSolid::GetName
G4String GetName() const
G4VSolid::ClipBetweenSections
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:376
G4VSolid::DumpInfo
void DumpInfo() const
G4VSolid::kCarTolerance
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4VSolid::ClipCrossSection
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:345
G4VoxelLimits::GetMinExtent
G4double GetMinExtent(const EAxis pAxis) const
G4VoxelLimits::IsYLimited
G4bool IsYLimited() const
G4VoxelLimits::GetMinZExtent
G4double GetMinZExtent() const
G4VoxelLimits::IsXLimited
G4bool IsXLimited() const
G4VoxelLimits::GetMaxExtent
G4double GetMaxExtent(const EAxis pAxis) const
G4VoxelLimits::GetMaxYExtent
G4double GetMaxYExtent() const
G4VoxelLimits::GetMaxZExtent
G4double GetMaxZExtent() const
G4VoxelLimits::GetMinYExtent
G4double GetMinYExtent() const
G4VoxelLimits::GetMinXExtent
G4double GetMinXExtent() const
G4VoxelLimits::IsZLimited
G4bool IsZLimited() const
G4VoxelLimits::GetMaxXExtent
G4double GetMaxXExtent() const
EAxis
EAxis
Definition: geomdefs.hh:54
kYAxis
@ kYAxis
Definition: geomdefs.hh:54
kXAxis
@ kXAxis
Definition: geomdefs.hh:54
kZAxis
@ kZAxis
Definition: geomdefs.hh:54
EInside
EInside
Definition: geomdefs.hh:58
kInside
@ kInside
Definition: geomdefs.hh:58
kOutside
@ kOutside
Definition: geomdefs.hh:58
kSurface
@ kSurface
Definition: geomdefs.hh:58
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
kMaxMeshSections
const G4int kMaxMeshSections
Definition: meshdefs.hh:46
kMinMeshSections
const G4int kMinMeshSections
Definition: meshdefs.hh:45
kMeshAngleDefault
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42
CLHEP
Definition: DoubConv.h:17
sqr
T sqr(const T &x)
Definition: templates.hh:145