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