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