source: trunk/source/geometry/solids/CSG/src/G4Cons.cc@ 1347

Last change on this file since 1347 was 1337, checked in by garnier, 15 years ago

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