| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // 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
55 using namespace CLHEP;
56
57 ////////////////////////////////////////////////////////////////////////
58 //
59 // Private enums: Not for external use
60
61 namespace
62 {
63 // used by DistanceToOut()
64 //
65 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
66
67 // used by ApproxSurfaceNormal()
68 //
69 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
70 }
71
72 //////////////////////////////////////////////////////////////////////////
73 //
74 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
75 // - note if pDPhi>2PI then reset to 2PI
76
77 G4Cons::G4Cons( const G4String& pName,
78 G4double pRmin1, G4double pRmax1,
79 G4double pRmin2, G4double pRmax2,
80 G4double pDz,
81 G4double pSPhi, G4double pDPhi)
82 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
83 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
84 {
85 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
86 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
87
88 halfCarTolerance=kCarTolerance*0.5;
89 halfRadTolerance=kRadTolerance*0.5;
90 halfAngTolerance=kAngTolerance*0.5;
91
92 // Check z-len
93 //
94 if ( pDz < 0 )
95 {
96 std::ostringstream message;
97 message << "Invalid Z half-length for Solid: " << GetName() << G4endl
98 << " hZ = " << pDz;
99 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
100 FatalException, message);
101 }
102
103 // Check radii
104 //
105 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
106 {
107 std::ostringstream message;
108 message << "Invalid values of radii for Solid: " << GetName() << G4endl
109 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
110 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
111 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
112 FatalException, message) ;
113 }
114 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
115 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
116
117 // Check angles
118 //
119 CheckPhiAngles(pSPhi, pDPhi);
120 }
121
122 ///////////////////////////////////////////////////////////////////////
123 //
124 // Fake default constructor - sets only member data and allocates memory
125 // for usage restricted to object persistency.
126 //
127 G4Cons::G4Cons( __void__& a )
128 : G4CSGSolid(a)
129 {
130 }
131
132 ///////////////////////////////////////////////////////////////////////
133 //
134 // Destructor
135
136 G4Cons::~G4Cons() = default;
137
138 //////////////////////////////////////////////////////////////////////////
139 //
140 // Copy constructor
141
142 G4Cons::G4Cons(const G4Cons&) = default;
143
144 //////////////////////////////////////////////////////////////////////////
145 //
146 // Assignment operator
147
148 G4Cons& G4Cons::operator = (const G4Cons& rhs)
149 {
150 // Check assignment to self
151 //
152 if (this == &rhs) { return *this; }
153
154 // Copy base class data
155 //
156 G4CSGSolid::operator=(rhs);
157
158 // Copy data
159 //
160 kRadTolerance = rhs.kRadTolerance;
161 kAngTolerance = rhs.kAngTolerance;
162 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
163 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
164 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
165 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; cosHDPhi = rhs.cosHDPhi;
166 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
167 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
168 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
169 fPhiFullCone = rhs.fPhiFullCone;
170 halfCarTolerance = rhs.halfCarTolerance;
171 halfRadTolerance = rhs.halfRadTolerance;
172 halfAngTolerance = rhs.halfAngTolerance;
173
174 return *this;
175 }
176
177 /////////////////////////////////////////////////////////////////////
178 //
179 // Return whether point inside/outside/on surface
180
181 EInside G4Cons::Inside(const G4ThreeVector& p) const
182 {
183 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
184 EInside in;
185
186 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
187 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
188 else { in = kInside; }
189
190 r2 = p.x()*p.x() + p.y()*p.y() ;
191 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
192 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
193
194 // rh2 = rh*rh;
195
196 tolRMin = rl - halfRadTolerance;
197 if ( tolRMin < 0 ) { tolRMin = 0; }
198 tolRMax = rh + halfRadTolerance;
199
200 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
201
202 if (rl != 0.0) { tolRMin = rl + halfRadTolerance; }
203 else { tolRMin = 0.0; }
204 tolRMax = rh - halfRadTolerance;
205
206 if (in == kInside) // else it's kSurface already
207 {
208 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
209 }
210 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
211 {
212 pPhi = std::atan2(p.y(),p.x()) ;
213
214 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
215 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
216
217 if ( (pPhi < fSPhi - halfAngTolerance) ||
218 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
219
220 else if (in == kInside) // else it's kSurface anyway already
221 {
222 if ( (pPhi < fSPhi + halfAngTolerance) ||
223 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
224 }
225 }
226 else if ( !fPhiFullCone ) { in = kSurface; }
227
228 return in ;
229 }
230
231 /////////////////////////////////////////////////////////////////////////
232 //
233 // Dispatch to parameterisation for replication mechanism dimension
234 // computation & modification.
235
236 void G4Cons::ComputeDimensions( G4VPVParameterisation* p,
237 const G4int n,
238 const G4VPhysicalVolume* pRep )
239 {
240 p->ComputeDimensions(*this,n,pRep) ;
241 }
242
243 ///////////////////////////////////////////////////////////////////////
244 //
245 // Get bounding box
246
247 void G4Cons::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
248 {
249 G4double rmin = std::min(GetInnerRadiusMinusZ(),GetInnerRadiusPlusZ());
250 G4double rmax = std::max(GetOuterRadiusMinusZ(),GetOuterRadiusPlusZ());
251 G4double dz = GetZHalfLength();
252
253 // Find bounding box
254 //
255 if (GetDeltaPhiAngle() < twopi)
256 {
257 G4TwoVector vmin,vmax;
258 G4GeomTools::DiskExtent(rmin,rmax,
259 GetSinStartPhi(),GetCosStartPhi(),
260 GetSinEndPhi(),GetCosEndPhi(),
261 vmin,vmax);
262 pMin.set(vmin.x(),vmin.y(),-dz);
263 pMax.set(vmax.x(),vmax.y(), dz);
264 }
265 else
266 {
267 pMin.set(-rmax,-rmax,-dz);
268 pMax.set( rmax, rmax, dz);
269 }
270
271 // Check correctness of the bounding box
272 //
273 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
274 {
275 std::ostringstream message;
276 message << "Bad bounding box (min >= max) for solid: "
277 << GetName() << " !"
278 << "\npMin = " << pMin
279 << "\npMax = " << pMax;
280 G4Exception("G4Cons::BoundingLimits()", "GeomMgt0001",
281 JustWarning, message);
282 DumpInfo();
283 }
284 }
285
286 ///////////////////////////////////////////////////////////////////////
287 //
288 // Calculate extent under transform and specified limit
289
290 G4bool G4Cons::CalculateExtent( const EAxis pAxis,
291 const G4VoxelLimits& pVoxelLimit,
292 const G4AffineTransform& pTransform,
293 G4double& pMin,
294 G4double& pMax ) const
295 {
296 G4ThreeVector bmin, bmax;
297 G4bool exist;
298
299 // Get bounding box
300 BoundingLimits(bmin,bmax);
301
302 // Check bounding box
303 G4BoundingEnvelope bbox(bmin,bmax);
304 #ifdef G4BBOX_EXTENT
305 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
306 #endif
307 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
308 {
309 return exist = pMin < pMax;
310 }
311
312 // Get parameters of the solid
313 G4double rmin1 = GetInnerRadiusMinusZ();
314 G4double rmax1 = GetOuterRadiusMinusZ();
315 G4double rmin2 = GetInnerRadiusPlusZ();
316 G4double rmax2 = GetOuterRadiusPlusZ();
317 G4double dz = GetZHalfLength();
318 G4double dphi = GetDeltaPhiAngle();
319
320 // Find bounding envelope and calculate extent
321 //
322 const G4int NSTEPS = 24; // number of steps for whole circle
323 G4double astep = twopi/NSTEPS; // max angle for one step
324 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
325 G4double ang = dphi/ksteps;
326
327 G4double sinHalf = std::sin(0.5*ang);
328 G4double cosHalf = std::cos(0.5*ang);
329 G4double sinStep = 2.*sinHalf*cosHalf;
330 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
331 G4double rext1 = rmax1/cosHalf;
332 G4double rext2 = rmax2/cosHalf;
333
334 // bounding envelope for full cone without hole consists of two polygons,
335 // in other cases it is a sequence of quadrilaterals
336 if (rmin1 == 0 && rmin2 == 0 && dphi == twopi)
337 {
338 G4double sinCur = sinHalf;
339 G4double cosCur = cosHalf;
340
341 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
342 for (G4int k=0; k<NSTEPS; ++k)
343 {
344 baseA[k].set(rext1*cosCur,rext1*sinCur,-dz);
345 baseB[k].set(rext2*cosCur,rext2*sinCur, dz);
346
347 G4double sinTmp = sinCur;
348 sinCur = sinCur*cosStep + cosCur*sinStep;
349 cosCur = cosCur*cosStep - sinTmp*sinStep;
350 }
351 std::vector<const G4ThreeVectorList *> polygons(2);
352 polygons[0] = &baseA;
353 polygons[1] = &baseB;
354 G4BoundingEnvelope benv(bmin,bmax,polygons);
355 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
356 }
357 else
358 {
359 G4double sinStart = GetSinStartPhi();
360 G4double cosStart = GetCosStartPhi();
361 G4double sinEnd = GetSinEndPhi();
362 G4double cosEnd = GetCosEndPhi();
363 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
364 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
365
366 // set quadrilaterals
367 G4ThreeVectorList pols[NSTEPS+2];
368 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
369 pols[0][0].set(rmin2*cosStart,rmin2*sinStart, dz);
370 pols[0][1].set(rmin1*cosStart,rmin1*sinStart,-dz);
371 pols[0][2].set(rmax1*cosStart,rmax1*sinStart,-dz);
372 pols[0][3].set(rmax2*cosStart,rmax2*sinStart, dz);
373 for (G4int k=1; k<ksteps+1; ++k)
374 {
375 pols[k][0].set(rmin2*cosCur,rmin2*sinCur, dz);
376 pols[k][1].set(rmin1*cosCur,rmin1*sinCur,-dz);
377 pols[k][2].set(rext1*cosCur,rext1*sinCur,-dz);
378 pols[k][3].set(rext2*cosCur,rext2*sinCur, dz);
379
380 G4double sinTmp = sinCur;
381 sinCur = sinCur*cosStep + cosCur*sinStep;
382 cosCur = cosCur*cosStep - sinTmp*sinStep;
383 }
384 pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*sinEnd, dz);
385 pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*sinEnd,-dz);
386 pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*sinEnd,-dz);
387 pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*sinEnd, dz);
388
389 // set envelope and calculate extent
390 std::vector<const G4ThreeVectorList *> polygons;
391 polygons.resize(ksteps+2);
392 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
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
405 G4ThreeVector 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
509 G4ThreeVector 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
647 G4double 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 else
847 {
848 xi = p.x() + sd*v.x() ;
849 yi = p.y() + sd*v.y() ;
850 ri = rMaxAv + zi*tanRMax ;
851 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
852
853 if ( cosPsi >= cosHDPhiIT ) { return sd; }
854 }
855 }
856 } // end if (sd>0)
857 }
858 }
859 else
860 {
861 // Inside outer cone
862 // check not inside, and heading through G4Cons (-> 0 to in)
863
864 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
865 (rin + halfRadTolerance*secRMin) )
866 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
867 {
868 // Inside cones, delta r -ve, inside z extent
869 // Point is on the Surface => check Direction using Normal.dot(v)
870
871 xi = p.x() ;
872 yi = p.y() ;
873 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
874 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
875 if ( !fPhiFullCone )
876 {
877 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
878 if ( cosPsi >= cosHDPhiIT )
879 {
880 if ( Normal.dot(v) <= 0 ) { return 0.0; }
881 }
882 }
883 else
884 {
885 if ( Normal.dot(v) <= 0 ) { return 0.0; }
886 }
887 }
888 }
889 }
890 else // Single root case
891 {
892 if ( std::fabs(nt2) > kRadTolerance )
893 {
894 sd = -0.5*nt3/nt2 ;
895
896 if ( sd < 0 ) { return kInfinity; } // travel away
897 else // sd >= 0, If 'forwards'. Check z intersection
898 {
899 zi = p.z() + sd*v.z() ;
900
901 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
902 {
903 // Z ok. Check phi intersection if reqd
904
905 if ( fPhiFullCone ) { return sd; }
906 else
907 {
908 xi = p.x() + sd*v.x() ;
909 yi = p.y() + sd*v.y() ;
910 ri = rMaxAv + zi*tanRMax ;
911 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
912
913 if (cosPsi >= cosHDPhiIT) { return sd; }
914 }
915 }
916 }
917 }
918 else // travel || cone surface from its origin
919 {
920 sd = kInfinity ;
921 }
922 }
923
924 // Inner Cone Intersection
925 // o Space is divided into 3 areas:
926 // 1) Radius greater than real inner cone & imaginary cone & outside
927 // tolerance
928 // 2) Radius less than inner or imaginary cone & outside tolarance
929 // 3) Within tolerance of real or imaginary cones
930 // - Extra checks needed for 3's intersections
931 // => lots of duplicated code
932
933 if (rMinAv != 0.0)
934 {
935 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
936 nt2 = t2 - tanRMin*v.z()*rin ;
937 nt3 = t3 - rin*rin ;
938
939 if ( nt1 != 0.0 )
940 {
941 if ( nt3 > rin*kRadTolerance*secRMin )
942 {
943 // At radius greater than real & imaginary cones
944 // -> 2nd root, with zi check
945
946 b = nt2/nt1 ;
947 c = nt3/nt1 ;
948 d = b*b-c ;
949 if (d >= 0) // > 0
950 {
951 if(b>0){sd = c/( -b-std::sqrt(d));}
952 else {sd = -b + std::sqrt(d) ;}
953
954 if ( sd >= 0 ) // > 0
955 {
956 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
957 { // 64 bits systems. Split long distance and recompute
958 G4double fTerm = sd-std::fmod(sd,dRmax);
959 sd = fTerm + DistanceToIn(p+fTerm*v,v);
960 }
961 zi = p.z() + sd*v.z() ;
962
963 if ( std::fabs(zi) <= tolODz )
964 {
965 if ( !fPhiFullCone )
966 {
967 xi = p.x() + sd*v.x() ;
968 yi = p.y() + sd*v.y() ;
969 ri = rMinAv + zi*tanRMin ;
970 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
971
972 if (cosPsi >= cosHDPhiIT)
973 {
974 if ( sd > halfRadTolerance ) { snxt=sd; }
975 else
976 {
977 // Calculate a normal vector in order to check Direction
978
979 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
980 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
981 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
982 }
983 }
984 }
985 else
986 {
987 if ( sd > halfRadTolerance ) { return sd; }
988 else
989 {
990 // Calculate a normal vector in order to check Direction
991
992 xi = p.x() + sd*v.x() ;
993 yi = p.y() + sd*v.y() ;
994 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
995 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
996 if ( Normal.dot(v) <= 0 ) { return sd; }
997 }
998 }
999 }
1000 }
1001 }
1002 }
1003 else if ( nt3 < -rin*kRadTolerance*secRMin )
1004 {
1005 // Within radius of inner cone (real or imaginary)
1006 // -> Try 2nd root, with checking intersection is with real cone
1007 // -> If check fails, try 1st root, also checking intersection is
1008 // on real cone
1009
1010 b = nt2/nt1 ;
1011 c = nt3/nt1 ;
1012 d = b*b - c ;
1013
1014 if ( d >= 0 ) // > 0
1015 {
1016 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1017 else { sd = -b + std::sqrt(d); }
1018 zi = p.z() + sd*v.z() ;
1019 ri = rMinAv + zi*tanRMin ;
1020
1021 if ( ri > 0 )
1022 {
1023 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1024 {
1025 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1026 { // seen on 64 bits systems. Split and recompute
1027 G4double fTerm = sd-std::fmod(sd,dRmax);
1028 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1029 }
1030 if ( !fPhiFullCone )
1031 {
1032 xi = p.x() + sd*v.x() ;
1033 yi = p.y() + sd*v.y() ;
1034 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1035
1036 if (cosPsi >= cosHDPhiOT)
1037 {
1038 if ( sd > halfRadTolerance ) { snxt=sd; }
1039 else
1040 {
1041 // Calculate a normal vector in order to check Direction
1042
1043 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1044 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1045 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1046 }
1047 }
1048 }
1049 else
1050 {
1051 if( sd > halfRadTolerance ) { return sd; }
1052 else
1053 {
1054 // Calculate a normal vector in order to check Direction
1055
1056 xi = p.x() + sd*v.x() ;
1057 yi = p.y() + sd*v.y() ;
1058 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1059 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1060 if ( Normal.dot(v) <= 0 ) { return sd; }
1061 }
1062 }
1063 }
1064 }
1065 else
1066 {
1067 if (b>0) { sd = -b - std::sqrt(d); }
1068 else { sd = c/(-b+std::sqrt(d)); }
1069 zi = p.z() + sd*v.z() ;
1070 ri = rMinAv + zi*tanRMin ;
1071
1072 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1073 {
1074 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1075 { // seen on 64 bits systems. Split and recompute
1076 G4double fTerm = sd-std::fmod(sd,dRmax);
1077 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1078 }
1079 if ( !fPhiFullCone )
1080 {
1081 xi = p.x() + sd*v.x() ;
1082 yi = p.y() + sd*v.y() ;
1083 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1084
1085 if (cosPsi >= cosHDPhiIT)
1086 {
1087 if ( sd > halfRadTolerance ) { snxt=sd; }
1088 else
1089 {
1090 // Calculate a normal vector in order to check Direction
1091
1092 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1093 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1094 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1095 }
1096 }
1097 }
1098 else
1099 {
1100 if ( sd > halfRadTolerance ) { return sd; }
1101 else
1102 {
1103 // Calculate a normal vector in order to check Direction
1104
1105 xi = p.x() + sd*v.x() ;
1106 yi = p.y() + sd*v.y() ;
1107 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1108 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1109 if ( Normal.dot(v) <= 0 ) { return sd; }
1110 }
1111 }
1112 }
1113 }
1114 }
1115 }
1116 else
1117 {
1118 // Within kRadTol*0.5 of inner cone (real OR imaginary)
1119 // ----> Check not travelling through (=>0 to in)
1120 // ----> if not:
1121 // -2nd root with validity check
1122
1123 if ( std::fabs(p.z()) <= tolODz )
1124 {
1125 if ( nt2 > 0 )
1126 {
1127 // Inside inner real cone, heading outwards, inside z range
1128
1129 if ( !fPhiFullCone )
1130 {
1131 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
1132
1133 if (cosPsi >= cosHDPhiIT) { return 0.0; }
1134 }
1135 else { return 0.0; }
1136 }
1137 else
1138 {
1139 // Within z extent, but not travelling through
1140 // -> 2nd root or kInfinity if 1st root on imaginary cone
1141
1142 b = nt2/nt1 ;
1143 c = nt3/nt1 ;
1144 d = b*b - c ;
1145
1146 if ( d >= 0 ) // > 0
1147 {
1148 if (b>0) { sd = -b - std::sqrt(d); }
1149 else { sd = c/(-b+std::sqrt(d)); }
1150 zi = p.z() + sd*v.z() ;
1151 ri = rMinAv + zi*tanRMin ;
1152
1153 if ( ri > 0 ) // 2nd root
1154 {
1155 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1156 else { sd = -b + std::sqrt(d); }
1157
1158 zi = p.z() + sd*v.z() ;
1159
1160 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1161 {
1162 if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1163 { // seen on 64 bits systems. Split and recompute
1164 G4double fTerm = sd-std::fmod(sd,dRmax);
1165 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1166 }
1167 if ( !fPhiFullCone )
1168 {
1169 xi = p.x() + sd*v.x() ;
1170 yi = p.y() + sd*v.y() ;
1171 ri = rMinAv + zi*tanRMin ;
1172 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1173
1174 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1175 }
1176 else { return sd; }
1177 }
1178 }
1179 else { return kInfinity; }
1180 }
1181 }
1182 }
1183 else // 2nd root
1184 {
1185 b = nt2/nt1 ;
1186 c = nt3/nt1 ;
1187 d = b*b - c ;
1188
1189 if ( d > 0 )
1190 {
1191 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1192 else { sd = -b + std::sqrt(d) ; }
1193 zi = p.z() + sd*v.z() ;
1194
1195 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1196 {
1197 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1198 { // seen on 64 bits systems. Split and recompute
1199 G4double fTerm = sd-std::fmod(sd,dRmax);
1200 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1201 }
1202 if ( !fPhiFullCone )
1203 {
1204 xi = p.x() + sd*v.x();
1205 yi = p.y() + sd*v.y();
1206 ri = rMinAv + zi*tanRMin ;
1207 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1208
1209 if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1210 }
1211 else { return sd; }
1212 }
1213 }
1214 }
1215 }
1216 }
1217 }
1218
1219 // Phi segment intersection
1220 //
1221 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1222 //
1223 // o NOTE: Large duplication of code between sphi & ephi checks
1224 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1225 // intersection check <=0 -> >=0
1226 // -> Should use some form of loop Construct
1227
1228 if ( !fPhiFullCone )
1229 {
1230 // First phi surface (starting phi)
1231
1232 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1233
1234 if ( Comp < 0 ) // Component in outwards normal dirn
1235 {
1236 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1237
1238 if (Dist < halfCarTolerance)
1239 {
1240 sd = Dist/Comp ;
1241
1242 if ( sd < snxt )
1243 {
1244 if ( sd < 0 ) { sd = 0.0; }
1245
1246 zi = p.z() + sd*v.z() ;
1247
1248 if ( std::fabs(zi) <= tolODz )
1249 {
1250 xi = p.x() + sd*v.x() ;
1251 yi = p.y() + sd*v.y() ;
1252 rhoi2 = xi*xi + yi*yi ;
1253 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1254 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1255
1256 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1257 {
1258 // z and r intersections good - check intersecting with
1259 // correct half-plane
1260
1261 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1262 }
1263 }
1264 }
1265 }
1266 }
1267
1268 // Second phi surface (Ending phi)
1269
1270 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1271
1272 if ( Comp < 0 ) // Component in outwards normal dirn
1273 {
1274 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1275 if (Dist < halfCarTolerance)
1276 {
1277 sd = Dist/Comp ;
1278
1279 if ( sd < snxt )
1280 {
1281 if ( sd < 0 ) { sd = 0.0; }
1282
1283 zi = p.z() + sd*v.z() ;
1284
1285 if (std::fabs(zi) <= tolODz)
1286 {
1287 xi = p.x() + sd*v.x() ;
1288 yi = p.y() + sd*v.y() ;
1289 rhoi2 = xi*xi + yi*yi ;
1290 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1291 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1292
1293 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1294 {
1295 // z and r intersections good - check intersecting with
1296 // correct half-plane
1297
1298 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1299 }
1300 }
1301 }
1302 }
1303 }
1304 }
1305 if (snxt < halfCarTolerance) { snxt = 0.; }
1306
1307 return snxt ;
1308 }
1309
1310 //////////////////////////////////////////////////////////////////////////////
1311 //
1312 // Calculate distance (<= actual) to closest surface of shape from outside
1313 // - Calculate distance to z, radial planes
1314 // - Only to phi planes if outside phi extent
1315 // - Return 0 if point inside
1316
1317 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
1318 {
1319 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1320 G4double tanRMin, secRMin, pRMin ;
1321 G4double tanRMax, secRMax, pRMax ;
1322
1323 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1324 safeZ = std::fabs(p.z()) - fDz ;
1325
1326 if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) )
1327 {
1328 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1329 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1330 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1331 safeR1 = (pRMin - rho)/secRMin ;
1332
1333 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1334 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1335 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1336 safeR2 = (rho - pRMax)/secRMax ;
1337
1338 if ( safeR1 > safeR2) { safe = safeR1; }
1339 else { safe = safeR2; }
1340 }
1341 else
1342 {
1343 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1344 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1345 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1346 safe = (rho - pRMax)/secRMax ;
1347 }
1348 if ( safeZ > safe ) { safe = safeZ; }
1349
1350 if ( !fPhiFullCone && (rho != 0.0) )
1351 {
1352 // Psi=angle from central phi to point
1353
1354 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1355
1356 if ( cosPsi < cosHDPhi ) // Point lies outside phi range
1357 {
1358 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1359 {
1360 safePhi = std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
1361 }
1362 else
1363 {
1364 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1365 }
1366 if ( safePhi > safe ) { safe = safePhi; }
1367 }
1368 }
1369 if ( safe < 0.0 ) { safe = 0.0; }
1370
1371 return safe ;
1372 }
1373
1374 ///////////////////////////////////////////////////////////////
1375 //
1376 // Calculate distance to surface of shape from 'inside', allowing for tolerance
1377 // - Only Calc rmax intersection if no valid rmin intersection
1378
1379 G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
1380 const G4ThreeVector& v,
1381 const G4bool calcNorm,
1382 G4bool* validNorm,
1383 G4ThreeVector* n) const
1384 {
1385 ESide side = kNull, sider = kNull, sidephi = kNull;
1386
1387 G4double snxt,srd,sphi,pdist ;
1388
1389 G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone
1390 G4double tanRMin, secRMin, rMinAv ; // Data for inner cone
1391
1392 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1393 G4double b, c, d, sr2, sr3 ;
1394
1395 // Vars for intersection within tolerance
1396
1397 ESide sidetol = kNull ;
1398 G4double slentol = kInfinity ;
1399
1400 // Vars for phi intersection:
1401
1402 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1403 G4double zi, ri, deltaRoi2 ;
1404
1405 // Z plane intersection
1406
1407 if ( v.z() > 0.0 )
1408 {
1409 pdist = fDz - p.z() ;
1410
1411 if (pdist > halfCarTolerance)
1412 {
1413 snxt = pdist/v.z() ;
1414 side = kPZ ;
1415 }
1416 else
1417 {
1418 if (calcNorm)
1419 {
1420 *n = G4ThreeVector(0,0,1) ;
1421 *validNorm = true ;
1422 }
1423 return snxt = 0.0;
1424 }
1425 }
1426 else if ( v.z() < 0.0 )
1427 {
1428 pdist = fDz + p.z() ;
1429
1430 if ( pdist > halfCarTolerance)
1431 {
1432 snxt = -pdist/v.z() ;
1433 side = kMZ ;
1434 }
1435 else
1436 {
1437 if ( calcNorm )
1438 {
1439 *n = G4ThreeVector(0,0,-1) ;
1440 *validNorm = true ;
1441 }
1442 return snxt = 0.0 ;
1443 }
1444 }
1445 else // Travel perpendicular to z axis
1446 {
1447 snxt = kInfinity ;
1448 side = kNull ;
1449 }
1450
1451 // Radial Intersections
1452 //
1453 // Intersection with outer cone (possible return) and
1454 // inner cone (must also check phi)
1455 //
1456 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1457 //
1458 // Intersects with x^2+y^2=(a*z+b)^2
1459 //
1460 // where a=tanRMax or tanRMin
1461 // b=rMaxAv or rMinAv
1462 //
1463 // (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 ;
1464 // t1 t2 t3
1465 //
1466 // \--------u-------/ \-----------v----------/ \---------w--------/
1467
1468 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1469 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1470 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
1471
1472
1473 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1474 t2 = p.x()*v.x() + p.y()*v.y() ;
1475 t3 = p.x()*p.x() + p.y()*p.y() ;
1476 rout = tanRMax*p.z() + rMaxAv ;
1477
1478 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1479 nt2 = t2 - tanRMax*v.z()*rout ;
1480 nt3 = t3 - rout*rout ;
1481
1482 if (v.z() > 0.0)
1483 {
1484 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1485 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1486 }
1487 else if (v.z() < 0.0)
1488 {
1489 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1490 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1491 }
1492 else
1493 {
1494 deltaRoi2 = 1.0;
1495 }
1496
1497 if ( (nt1 != 0.0) && (deltaRoi2 > 0.0) )
1498 {
1499 // Equation quadratic => 2 roots : second root must be leaving
1500
1501 b = nt2/nt1 ;
1502 c = nt3/nt1 ;
1503 d = b*b - c ;
1504
1505 if ( d >= 0 )
1506 {
1507 // Check if on outer cone & heading outwards
1508 // NOTE: Should use rho-rout>-kRadTolerance*0.5
1509
1510 if (nt3 > -halfRadTolerance && nt2 >= 0 )
1511 {
1512 if (calcNorm)
1513 {
1514 risec = std::sqrt(t3)*secRMax ;
1515 *validNorm = true ;
1516 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1517 }
1518 return snxt=0 ;
1519 }
1520 else
1521 {
1522 sider = kRMax ;
1523 if (b>0) { srd = -b - std::sqrt(d); }
1524 else { srd = c/(-b+std::sqrt(d)) ; }
1525
1526 zi = p.z() + srd*v.z() ;
1527 ri = tanRMax*zi + rMaxAv ;
1528
1529 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1530 {
1531 // An intersection within the tolerance
1532 // we will Store it in case it is good -
1533 //
1534 slentol = srd ;
1535 sidetol = kRMax ;
1536 }
1537 if ( (ri < 0) || (srd < halfRadTolerance) )
1538 {
1539 // Safety: if both roots -ve ensure that srd cannot `win'
1540 // distance to out
1541
1542 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1543 else { sr2 = -b + std::sqrt(d); }
1544 zi = p.z() + sr2*v.z() ;
1545 ri = tanRMax*zi + rMaxAv ;
1546
1547 if ((ri >= 0) && (sr2 > halfRadTolerance))
1548 {
1549 srd = sr2;
1550 }
1551 else
1552 {
1553 srd = kInfinity ;
1554
1555 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
1556 {
1557 // An intersection within the tolerance.
1558 // Storing it in case it is good.
1559
1560 slentol = sr2 ;
1561 sidetol = kRMax ;
1562 }
1563 }
1564 }
1565 }
1566 }
1567 else
1568 {
1569 // No intersection with outer cone & not parallel
1570 // -> already outside, no intersection
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 }
1581 else if ( (nt2 != 0.0) && (deltaRoi2 > 0.0) )
1582 {
1583 // Linear case (only one intersection) => point outside outer cone
1584
1585 if ( calcNorm )
1586 {
1587 risec = std::sqrt(t3)*secRMax;
1588 *validNorm = true;
1589 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1590 }
1591 return snxt = 0.0 ;
1592 }
1593 else
1594 {
1595 // No intersection -> parallel to outer cone
1596 // => Z or inner cone intersection
1597
1598 srd = kInfinity ;
1599 }
1600
1601 // Check possible intersection within tolerance
1602
1603 if ( slentol <= halfCarTolerance )
1604 {
1605 // An intersection within the tolerance was found.
1606 // We must accept it only if the momentum points outwards.
1607 //
1608 // G4ThreeVector ptTol ; // The point of the intersection
1609 // ptTol= p + slentol*v ;
1610 // ri=tanRMax*zi+rMaxAv ;
1611 //
1612 // Calculate a normal vector, as below
1613
1614 xi = p.x() + slentol*v.x();
1615 yi = p.y() + slentol*v.y();
1616 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1617 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
1618
1619 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly
1620 {
1621 if ( calcNorm )
1622 {
1623 *n = Normal.unit() ;
1624 *validNorm = true ;
1625 }
1626 return snxt = 0.0 ;
1627 }
1628 else // On the surface, but not heading out so we ignore this intersection
1629 { // (as it is within tolerance).
1630 slentol = kInfinity ;
1631 }
1632 }
1633
1634 // Inner Cone intersection
1635
1636 if ( (fRmin1 != 0.0) || (fRmin2 != 0.0) )
1637 {
1638 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1639 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1640
1641 if ( nt1 != 0.0 )
1642 {
1643 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1644 rMinAv = (fRmin1 + fRmin2)*0.5 ;
1645 rin = tanRMin*p.z() + rMinAv ;
1646 nt2 = t2 - tanRMin*v.z()*rin ;
1647 nt3 = t3 - rin*rin ;
1648
1649 // Equation quadratic => 2 roots : first root must be leaving
1650
1651 b = nt2/nt1 ;
1652 c = nt3/nt1 ;
1653 d = b*b - c ;
1654
1655 if ( d >= 0.0 )
1656 {
1657 // NOTE: should be rho-rin<kRadTolerance*0.5,
1658 // but using squared versions for efficiency
1659
1660 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
1661 {
1662 if ( nt2 < 0.0 )
1663 {
1664 if (calcNorm) { *validNorm = false; }
1665 return snxt = 0.0;
1666 }
1667 }
1668 else
1669 {
1670 if (b>0) { sr2 = -b - std::sqrt(d); }
1671 else { sr2 = c/(-b+std::sqrt(d)); }
1672 zi = p.z() + sr2*v.z() ;
1673 ri = tanRMin*zi + rMinAv ;
1674
1675 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
1676 {
1677 // An intersection within the tolerance
1678 // storing it in case it is good.
1679
1680 slentol = sr2 ;
1681 sidetol = kRMax ;
1682 }
1683 if( (ri<0) || (sr2 < halfRadTolerance) )
1684 {
1685 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1686 else { sr3 = -b + std::sqrt(d) ; }
1687
1688 // Safety: if both roots -ve ensure that srd cannot `win'
1689 // distancetoout
1690
1691 if ( sr3 > halfRadTolerance )
1692 {
1693 if( sr3 < srd )
1694 {
1695 zi = p.z() + sr3*v.z() ;
1696 ri = tanRMin*zi + rMinAv ;
1697
1698 if ( ri >= 0.0 )
1699 {
1700 srd=sr3 ;
1701 sider=kRMin ;
1702 }
1703 }
1704 }
1705 else if ( sr3 > -halfRadTolerance )
1706 {
1707 // Intersection in tolerance. Store to check if it's good
1708
1709 slentol = sr3 ;
1710 sidetol = kRMin ;
1711 }
1712 }
1713 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
1714 {
1715 srd = sr2 ;
1716 sider = kRMin ;
1717 }
1718 else if (sr2 > -halfCarTolerance)
1719 {
1720 // Intersection in tolerance. Store to check if it's good
1721
1722 slentol = sr2 ;
1723 sidetol = kRMin ;
1724 }
1725 if( slentol <= halfCarTolerance )
1726 {
1727 // An intersection within the tolerance was found.
1728 // We must accept it only if the momentum points outwards.
1729
1730 G4ThreeVector Normal ;
1731
1732 // Calculate a normal vector, as below
1733
1734 xi = p.x() + slentol*v.x() ;
1735 yi = p.y() + slentol*v.y() ;
1736 if( sidetol==kRMax )
1737 {
1738 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1739 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1740 }
1741 else
1742 {
1743 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1744 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1745 }
1746 if( Normal.dot(v) > 0 )
1747 {
1748 // We will leave the cone immediately
1749
1750 if( calcNorm )
1751 {
1752 *n = Normal.unit() ;
1753 *validNorm = true ;
1754 }
1755 return snxt = 0.0 ;
1756 }
1757 else
1758 {
1759 // On the surface, but not heading out so we ignore this
1760 // intersection (as it is within tolerance).
1761
1762 slentol = kInfinity ;
1763 }
1764 }
1765 }
1766 }
1767 }
1768 }
1769
1770 // Linear case => point outside inner cone ---> outer cone intersect
1771 //
1772 // Phi Intersection
1773
1774 if ( !fPhiFullCone )
1775 {
1776 // add angle calculation with correction
1777 // of the difference in domain of atan2 and Sphi
1778
1779 vphi = std::atan2(v.y(),v.x()) ;
1780
1781 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1782 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1783
1784 if ( (p.x() != 0.0) || (p.y() != 0.0) ) // Check if on z axis (rho not needed later)
1785 {
1786 // pDist -ve when inside
1787
1788 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1789 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1790
1791 // Comp -ve when in direction of outwards normal
1792
1793 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1794 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1795
1796 sidephi = kNull ;
1797
1798 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1799 && (pDistE <= halfCarTolerance) ) )
1800 || ( (fDPhi > pi) && ((pDistS <= halfCarTolerance)
1801 || (pDistE <= halfCarTolerance) ) ) )
1802 {
1803 // Inside both phi *full* planes
1804 if ( compS < 0 )
1805 {
1806 sphi = pDistS/compS ;
1807 if (sphi >= -halfCarTolerance)
1808 {
1809 xi = p.x() + sphi*v.x() ;
1810 yi = p.y() + sphi*v.y() ;
1811
1812 // Check intersecting with correct half-plane
1813 // (if not -> no intersect)
1814 //
1815 if ( (std::fabs(xi)<=kCarTolerance)
1816 && (std::fabs(yi)<=kCarTolerance) )
1817 {
1818 sidephi= kSPhi;
1819 if ( ( fSPhi-halfAngTolerance <= vphi )
1820 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
1821 {
1822 sphi = kInfinity;
1823 }
1824 }
1825 else
1826 if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
1827 {
1828 sphi = kInfinity ;
1829 }
1830 else
1831 {
1832 sidephi = kSPhi ;
1833 if ( pDistS > -halfCarTolerance )
1834 {
1835 sphi = 0.0 ; // Leave by sphi immediately
1836 }
1837 }
1838 }
1839 else
1840 {
1841 sphi = kInfinity ;
1842 }
1843 }
1844 else
1845 {
1846 sphi = kInfinity ;
1847 }
1848
1849 if ( compE < 0 )
1850 {
1851 sphi2 = pDistE/compE ;
1852
1853 // Only check further if < starting phi intersection
1854 //
1855 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1856 {
1857 xi = p.x() + sphi2*v.x() ;
1858 yi = p.y() + sphi2*v.y() ;
1859
1860 // Check intersecting with correct half-plane
1861
1862 if ( (std::fabs(xi)<=kCarTolerance)
1863 && (std::fabs(yi)<=kCarTolerance) )
1864 {
1865 // Leaving via ending phi
1866
1867 if( (fSPhi-halfAngTolerance > vphi)
1868 || (fSPhi+fDPhi+halfAngTolerance < vphi) )
1869 {
1870 sidephi = kEPhi ;
1871 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1872 else { sphi = 0.0; }
1873 }
1874 }
1875 else // Check intersecting with correct half-plane
1876 if ( yi*cosCPhi-xi*sinCPhi >= 0 )
1877 {
1878 // Leaving via ending phi
1879
1880 sidephi = kEPhi ;
1881 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1882 else { sphi = 0.0; }
1883 }
1884 }
1885 }
1886 }
1887 else
1888 {
1889 sphi = kInfinity ;
1890 }
1891 }
1892 else
1893 {
1894 // On z axis + travel not || to z axis -> if phi of vector direction
1895 // within phi of shape, Step limited by rmax, else Step =0
1896
1897 if ( (fSPhi-halfAngTolerance <= vphi)
1898 && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
1899 {
1900 sphi = kInfinity ;
1901 }
1902 else
1903 {
1904 sidephi = kSPhi ; // arbitrary
1905 sphi = 0.0 ;
1906 }
1907 }
1908 if ( sphi < snxt ) // Order intersecttions
1909 {
1910 snxt = sphi ;
1911 side = sidephi ;
1912 }
1913 }
1914 if ( srd < snxt ) // Order intersections
1915 {
1916 snxt = srd ;
1917 side = sider ;
1918 }
1919 if (calcNorm)
1920 {
1921 switch(side)
1922 { // Note: returned vector not normalised
1923 case kRMax: // (divide by frmax for unit vector)
1924 xi = p.x() + snxt*v.x() ;
1925 yi = p.y() + snxt*v.y() ;
1926 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1927 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1928 *validNorm = true ;
1929 break ;
1930 case kRMin:
1931 *validNorm = false ; // Rmin is inconvex
1932 break ;
1933 case kSPhi:
1934 if ( fDPhi <= pi )
1935 {
1936 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0);
1937 *validNorm = true ;
1938 }
1939 else
1940 {
1941 *validNorm = false ;
1942 }
1943 break ;
1944 case kEPhi:
1945 if ( fDPhi <= pi )
1946 {
1947 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
1948 *validNorm = true ;
1949 }
1950 else
1951 {
1952 *validNorm = false ;
1953 }
1954 break ;
1955 case kPZ:
1956 *n = G4ThreeVector(0,0,1) ;
1957 *validNorm = true ;
1958 break ;
1959 case kMZ:
1960 *n = G4ThreeVector(0,0,-1) ;
1961 *validNorm = true ;
1962 break ;
1963 default:
1964 G4cout << G4endl ;
1965 DumpInfo();
1966 std::ostringstream message;
1967 G4long oldprc = message.precision(16) ;
1968 message << "Undefined side for valid surface normal to solid."
1969 << G4endl
1970 << "Position:" << G4endl << G4endl
1971 << "p.x() = " << p.x()/mm << " mm" << G4endl
1972 << "p.y() = " << p.y()/mm << " mm" << G4endl
1973 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1974 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
1975 << " mm" << G4endl << G4endl ;
1976 if( p.x() != 0. || p.y() != 0.)
1977 {
1978 message << "point phi = " << std::atan2(p.y(),p.x())/degree
1979 << " degrees" << G4endl << G4endl ;
1980 }
1981 message << "Direction:" << G4endl << G4endl
1982 << "v.x() = " << v.x() << G4endl
1983 << "v.y() = " << v.y() << G4endl
1984 << "v.z() = " << v.z() << G4endl<< G4endl
1985 << "Proposed distance :" << G4endl<< G4endl
1986 << "snxt = " << snxt/mm << " mm" << G4endl ;
1987 message.precision(oldprc) ;
1988 G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002",
1989 JustWarning, message) ;
1990 break ;
1991 }
1992 }
1993 if (snxt < halfCarTolerance) { snxt = 0.; }
1994
1995 return snxt ;
1996 }
1997
1998 //////////////////////////////////////////////////////////////////
1999 //
2000 // Calculate distance (<=actual) to closest surface of shape from inside
2001
2002 G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const
2003 {
2004 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2005 G4double tanRMin, secRMin, pRMin;
2006 G4double tanRMax, secRMax, pRMax;
2007
2008 #ifdef G4CSGDEBUG
2009 if( Inside(p) == kOutside )
2010 {
2011 G4int oldprc=G4cout.precision(16) ;
2012 G4cout << G4endl ;
2013 DumpInfo();
2014 G4cout << "Position:" << G4endl << G4endl ;
2015 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2016 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2017 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2018 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2019 << " mm" << G4endl << G4endl ;
2020 if( (p.x() != 0.) || (p.x() != 0.) )
2021 {
2022 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2023 << " degrees" << G4endl << G4endl ;
2024 }
2025 G4cout.precision(oldprc) ;
2026 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
2027 JustWarning, "Point p is outside !?" );
2028 }
2029 #endif
2030
2031 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2032 safeZ = fDz - std::fabs(p.z()) ;
2033
2034 if ((fRmin1 != 0.0) || (fRmin2 != 0.0))
2035 {
2036 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2037 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2038 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2039 safeR1 = (rho - pRMin)/secRMin ;
2040 }
2041 else
2042 {
2043 safeR1 = kInfinity ;
2044 }
2045
2046 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2047 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2048 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2049 safeR2 = (pRMax - rho)/secRMax ;
2050
2051 if (safeR1 < safeR2) { safe = safeR1; }
2052 else { safe = safeR2; }
2053 if (safeZ < safe) { safe = safeZ ; }
2054
2055 // Check if phi divided, Calc distances closest phi plane
2056
2057 if (!fPhiFullCone)
2058 {
2059 // Above/below central phi of G4Cons?
2060
2061 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2062 {
2063 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2064 }
2065 else
2066 {
2067 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2068 }
2069 if (safePhi < safe) { safe = safePhi; }
2070 }
2071 if ( safe < 0 ) { safe = 0; }
2072
2073 return safe ;
2074 }
2075
2076 //////////////////////////////////////////////////////////////////////////
2077 //
2078 // GetEntityType
2079
2080 G4GeometryType G4Cons::GetEntityType() const
2081 {
2082 return {"G4Cons"};
2083 }
2084
2085 //////////////////////////////////////////////////////////////////////////
2086 //
2087 // Make a clone of the object
2088 //
2089 G4VSolid* G4Cons::Clone() const
2090 {
2091 return new G4Cons(*this);
2092 }
2093
2094 //////////////////////////////////////////////////////////////////////////
2095 //
2096 // Stream object contents to an output stream
2097
2098 std::ostream& G4Cons::StreamInfo(std::ostream& os) const
2099 {
2100 G4long oldprc = os.precision(16);
2101 os << "-----------------------------------------------------------\n"
2102 << " *** Dump for solid - " << GetName() << " ***\n"
2103 << " ===================================================\n"
2104 << " Solid type: G4Cons\n"
2105 << " Parameters: \n"
2106 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2107 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2108 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2109 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2110 << " half length in Z : " << fDz/mm << " mm \n"
2111 << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2112 << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2113 << "-----------------------------------------------------------\n";
2114 os.precision(oldprc);
2115
2116 return os;
2117 }
2118
2119
2120
2121 /////////////////////////////////////////////////////////////////////////
2122 //
2123 // GetPointOnSurface
2124
2125 G4ThreeVector G4Cons::GetPointOnSurface() const
2126 {
2127 // declare working variables
2128 //
2129 G4double rone = (fRmax1-fRmax2)/(2.*fDz);
2130 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz);
2131 G4double qone = (fRmax1 == fRmax2) ? 0. : fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
2132 G4double qtwo = (fRmin1 == fRmin2) ? 0. : fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
2133
2134 G4double slin = std::hypot(fRmin1-fRmin2, 2.*fDz);
2135 G4double slout = std::hypot(fRmax1-fRmax2, 2.*fDz);
2136 G4double Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout; // outer surface
2137 G4double Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin; // inner surface
2138 G4double Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1); // base at -Dz
2139 G4double Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2); // base at +Dz
2140 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); // phi section
2141
2142 G4double phi = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi);
2143 G4double cosu = std::cos(phi);
2144 G4double sinu = std::sin(phi);
2145 G4double rRand1 = GetRadiusInRing(fRmin1, fRmax1);
2146 G4double rRand2 = GetRadiusInRing(fRmin2, fRmax2);
2147
2148 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
2149 G4double chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2150
2151 if( (chose >= 0.) && (chose < Aone) ) // outer surface
2152 {
2153 if(fRmax1 != fRmax2)
2154 {
2155 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2156 return { rone*cosu*(qone-zRand), rone*sinu*(qone-zRand), zRand };
2157 }
2158 else
2159 {
2160 return { fRmax1*cosu, fRmax2*sinu, G4RandFlat::shoot(-1.*fDz,fDz) };
2161 }
2162 }
2163 else if( (chose >= Aone) && (chose < Aone + Atwo) ) // inner surface
2164 {
2165 if(fRmin1 != fRmin2)
2166 {
2167 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2168 return { rtwo*cosu*(qtwo-zRand), rtwo*sinu*(qtwo-zRand), zRand };
2169 }
2170 else
2171 {
2172 return { fRmin1*cosu, fRmin2*sinu, G4RandFlat::shoot(-1.*fDz,fDz) };
2173 }
2174 }
2175 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) // base at -Dz
2176 {
2177 return {rRand1*cosu, rRand1*sinu, -1*fDz};
2178 }
2179 else if( (chose >= Aone + Atwo + Athree)
2180 && (chose < Aone + Atwo + Athree + Afour) ) // base at +Dz
2181 {
2182 return { rRand2*cosu, rRand2*sinu, fDz };
2183 }
2184 else if( (chose >= Aone + Atwo + Athree + Afour) // SPhi section
2185 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2186 {
2187 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2188 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2189 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2190 return { rRand1*cosSPhi, rRand1*sinSPhi, zRand };
2191 }
2192 else // SPhi+DPhi section
2193 {
2194 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2195 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2196 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2197 return { rRand1*cosEPhi, rRand1*sinEPhi, zRand };
2198 }
2199 }
2200
2201 //////////////////////////////////////////////////////////////////////////
2202 //
2203 // Methods for visualisation
2204
2205 void G4Cons::DescribeYourselfTo (G4VGraphicsScene& scene) const
2206 {
2207 scene.AddSolid (*this);
2208 }
2209
2210 G4Polyhedron* G4Cons::CreatePolyhedron () const
2211 {
2212 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
2213 }
2214
2215 #endif
2216