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

Last change on this file since 833 was 831, checked in by garnier, 17 years ago

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