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