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

Last change on this file since 1089 was 1058, checked in by garnier, 17 years ago

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