| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Implementation of G4EllipticalCone class
27 //
28 // This code implements an Elliptical Cone giv
29 // equation:
30 // x^2/a^2 + y^2/b^2 = (z-h)^2
31 // and specified by the parameters (a,b,h) and
32 // xy plane above z = 0.
33 //
34 // Author: Dionysios Anninos
35 // Revised: Evgueni Tcherniaev
36 // -------------------------------------------
37
38 #if !(defined(G4GEOM_USE_UELLIPTICALCONE) && d
39
40 #include "globals.hh"
41
42 #include "G4EllipticalCone.hh"
43
44 #include "G4RandomTools.hh"
45 #include "G4GeomTools.hh"
46 #include "G4ClippablePolygon.hh"
47 #include "G4VoxelLimits.hh"
48 #include "G4AffineTransform.hh"
49 #include "G4BoundingEnvelope.hh"
50 #include "G4GeometryTolerance.hh"
51
52 #include "meshdefs.hh"
53
54 #include "Randomize.hh"
55
56 #include "G4VGraphicsScene.hh"
57 #include "G4VisExtent.hh"
58
59 #include "G4AutoLock.hh"
60
61 namespace
62 {
63 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE
64 }
65
66 using namespace CLHEP;
67
68 //////////////////////////////////////////////
69 //
70 // Constructor - check parameters
71
72 G4EllipticalCone::G4EllipticalCone(const G4Str
73 G4dou
74 G4dou
75 G4dou
76 G4dou
77 : G4VSolid(pName), zTopCut(0.)
78 {
79 halfCarTol = 0.5*kCarTolerance;
80
81 // Check Semi-Axis & Z-cut
82 //
83 if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.
84 {
85 std::ostringstream message;
86 message << "Invalid semi-axis or height f
87 << "\n X semi-axis, Y semi-axis
88 << pxSemiAxis << ", " << pySemiAx
89 G4Exception("G4EllipticalCone::G4Elliptic
90 FatalErrorInArgument, message
91 }
92
93 if ( pzTopCut <= 0 )
94 {
95 std::ostringstream message;
96 message << "Invalid z-coordinate for cutt
97 << "\n Z top cut = " << pzTopCu
98 G4Exception("G4EllipticalCone::G4Elliptic
99 FatalErrorInArgument, message
100 }
101
102 SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax )
103 SetZCut(pzTopCut);
104 }
105
106 //////////////////////////////////////////////
107 //
108 // Fake default constructor - sets only member
109 // for usage restri
110
111 G4EllipticalCone::G4EllipticalCone( __void__&
112 : G4VSolid(a), halfCarTol(0.),
113 xSemiAxis(0.), ySemiAxis(0.), zheight(0.),
114 cosAxisMin(0.), invXX(0.), invYY(0.)
115 {
116 }
117
118 //////////////////////////////////////////////
119 //
120 // Destructor
121
122 G4EllipticalCone::~G4EllipticalCone()
123 {
124 delete fpPolyhedron; fpPolyhedron = nullptr;
125 }
126
127 //////////////////////////////////////////////
128 //
129 // Copy constructor
130
131 G4EllipticalCone::G4EllipticalCone(const G4Ell
132 : G4VSolid(rhs), halfCarTol(rhs.halfCarTol),
133 fCubicVolume(rhs.fCubicVolume), fSurfaceAr
134 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.yS
135 zheight(rhs.zheight), zTopCut(rhs.zTopCut)
136 cosAxisMin(rhs.cosAxisMin), invXX(rhs.invX
137 {
138 }
139
140 //////////////////////////////////////////////
141 //
142 // Assignment operator
143
144 G4EllipticalCone& G4EllipticalCone::operator =
145 {
146 // Check assignment to self
147 //
148 if (this == &rhs) { return *this; }
149
150 // Copy base class data
151 //
152 G4VSolid::operator=(rhs);
153
154 // Copy data
155 //
156 halfCarTol = rhs.halfCarTol;
157 fCubicVolume = rhs.fCubicVolume; fSurfaceAr
158 xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.
159 zheight = rhs.zheight; zTopCut = rhs.zTopCu
160 cosAxisMin = rhs.cosAxisMin; invXX = rhs.in
161
162 fRebuildPolyhedron = false;
163 delete fpPolyhedron; fpPolyhedron = nullptr
164
165 return *this;
166 }
167
168 //////////////////////////////////////////////
169 //
170 // Get bounding box
171
172 void G4EllipticalCone::BoundingLimits(G4ThreeV
173 G4ThreeV
174 {
175 G4double zcut = GetZTopCut();
176 G4double height = GetZMax();
177 G4double xmax = GetSemiAxisX()*(height+zcu
178 G4double ymax = GetSemiAxisY()*(height+zcu
179 pMin.set(-xmax,-ymax,-zcut);
180 pMax.set( xmax, ymax, zcut);
181
182 // Check correctness of the bounding box
183 //
184 if (pMin.x() >= pMax.x() || pMin.y() >= pMax
185 {
186 std::ostringstream message;
187 message << "Bad bounding box (min >= max)
188 << GetName() << " !"
189 << "\npMin = " << pMin
190 << "\npMax = " << pMax;
191 G4Exception("G4EllipticalCone::BoundingLim
192 JustWarning, message);
193 DumpInfo();
194 }
195 }
196
197 //////////////////////////////////////////////
198 //
199 // Calculate extent under transform and specif
200
201 G4bool
202 G4EllipticalCone::CalculateExtent(const EAxis
203 const G4Voxe
204 const G4Affi
205 G4doub
206 {
207 G4ThreeVector bmin,bmax;
208 G4bool exist;
209
210 // Check bounding box (bbox)
211 //
212 BoundingLimits(bmin,bmax);
213 G4BoundingEnvelope bbox(bmin,bmax);
214 #ifdef G4BBOX_EXTENT
215 return bbox.CalculateExtent(pAxis,pVoxelLimi
216 #endif
217 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox
218 {
219 return exist = pMin < pMax;
220 }
221
222 // Set bounding envelope (benv) and calculat
223 //
224 static const G4int NSTEPS = 48; // number of
225 static const G4double ang = twopi/NSTEPS;
226 static const G4double sinHalf = std::sin(0.5
227 static const G4double cosHalf = std::cos(0.5
228 static const G4double sinStep = 2.*sinHalf*c
229 static const G4double cosStep = 1. - 2.*sinH
230 G4double zcut = bmax.z();
231 G4double height = GetZMax();
232 G4double sxmin = GetSemiAxisX()*(height-zcu
233 G4double symin = GetSemiAxisY()*(height-zcu
234 G4double sxmax = bmax.x()/cosHalf;
235 G4double symax = bmax.y()/cosHalf;
236
237 G4double sinCur = sinHalf;
238 G4double cosCur = cosHalf;
239 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS
240 for (G4int k=0; k<NSTEPS; ++k)
241 {
242 baseA[k].set(sxmax*cosCur,symax*sinCur,-zc
243 baseB[k].set(sxmin*cosCur,symin*sinCur, zc
244
245 G4double sinTmp = sinCur;
246 sinCur = sinCur*cosStep + cosCur*sinStep;
247 cosCur = cosCur*cosStep - sinTmp*sinStep;
248 }
249
250 std::vector<const G4ThreeVectorList *> polyg
251 polygons[0] = &baseA;
252 polygons[1] = &baseB;
253 G4BoundingEnvelope benv(bmin,bmax,polygons);
254 exist = benv.CalculateExtent(pAxis,pVoxelLim
255 return exist;
256 }
257
258 //////////////////////////////////////////////
259 //
260 // Determine where is point: inside, outside o
261
262 EInside G4EllipticalCone::Inside(const G4Three
263 {
264 G4double hp = std::sqrt(p.x()*p.x()*invXX +
265 G4double ds = (hp - zheight)*cosAxisMin;
266 G4double dz = std::abs(p.z()) - zTopCut;
267 G4double dist = std::max(ds,dz);
268
269 if (dist > halfCarTol) return kOutside;
270 return (dist > -halfCarTol) ? kSurface : kIn
271 }
272
273 //////////////////////////////////////////////
274 //
275 // Return unit normal at surface closest to p
276
277 G4ThreeVector G4EllipticalCone::SurfaceNormal(
278 {
279 G4ThreeVector norm(0,0,0);
280 G4int nsurf = 0; // number of surfaces wher
281
282 G4double hp = std::sqrt(p.x()*p.x()*invXX +
283 G4double ds = (hp - zheight)*cosAxisMin;
284 if (std::abs(ds) <= halfCarTol)
285 {
286 norm = G4ThreeVector(p.x()*invXX, p.y()*in
287 G4double mag = norm.mag();
288 if (mag == 0) return {0,0,1}; // apex
289 norm *= (1/mag);
290 ++nsurf;
291 }
292 G4double dz = std::abs(p.z()) - zTopCut;
293 if (std::abs(dz) <= halfCarTol)
294 {
295 norm += G4ThreeVector(0., 0.,(p.z() < 0) ?
296 ++nsurf;
297 }
298
299 if (nsurf == 1) return norm;
300 else if (nsurf > 1) return norm.unit(); //
301 else
302 {
303 // Point is not on the surface
304 //
305 #ifdef G4CSGDEBUG
306 std::ostringstream message;
307 G4long oldprc = message.precision(16);
308 message << "Point p is not on surface (!?)
309 << GetName() << G4endl;
310 message << "Position:\n";
311 message << " p.x() = " << p.x()/mm << "
312 message << " p.y() = " << p.y()/mm << "
313 message << " p.z() = " << p.z()/mm << "
314 G4cout.precision(oldprc);
315 G4Exception("G4EllipticalCone::SurfaceNorm
316 JustWarning, message );
317 DumpInfo();
318 #endif
319 return ApproxSurfaceNormal(p);
320 }
321 }
322
323 //////////////////////////////////////////////
324 //
325 // Find surface nearest to point and return co
326 // The algorithm is similar to the algorithm u
327 // This method normally should not be called.
328
329 G4ThreeVector
330 G4EllipticalCone::ApproxSurfaceNormal(const G4
331 {
332 G4double hp = std::sqrt(p.x()*p.x()*invXX +
333 G4double ds = (hp - zheight)*cosAxisMin;
334 G4double dz = std::abs(p.z()) - zTopCut;
335 if (ds > dz && std::abs(hp - p.z()) > halfCa
336 return G4ThreeVector(p.x()*invXX, p.y()*in
337 else
338 return { 0., 0., (G4double)((p.z() < 0) ?
339 }
340
341 //////////////////////////////////////////////
342 //
343 // Calculate distance to shape from outside, a
344 // return kInfinity if no intersection, or int
345
346 G4double G4EllipticalCone::DistanceToIn( const
347 const
348 {
349 G4double distMin = kInfinity;
350
351 // code from EllipticalTube
352
353 G4double sigz = p.z()+zTopCut;
354
355 //
356 // Check z = -dz planer surface
357 //
358
359 if (sigz < halfCarTol)
360 {
361 //
362 // We are "behind" the shape in z, and so
363 // potentially hit the rear face. Correct
364 //
365 if (v.z() <= 0)
366 {
367 //
368 // As long as we are far enough away, we
369 // can't intersect
370 //
371 if (sigz < 0) return kInfinity;
372
373 //
374 // Otherwise, we don't intersect unless
375 // on the surface of the ellipse
376 //
377
378 if ( sqr(p.x()/( xSemiAxis - halfCarTol
379 + sqr(p.y()/( ySemiAxis - halfCarTol
380 return kInfinity;
381
382 }
383 else
384 {
385 //
386 // How far?
387 //
388 G4double q = -sigz/v.z();
389
390 //
391 // Where does that place us?
392 //
393 G4double xi = p.x() + q*v.x(),
394 yi = p.y() + q*v.y();
395
396 //
397 // Is this on the surface (within ellips
398 //
399 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxi
400 {
401 //
402 // Yup. Return q, unless we are on the
403 //
404 return (sigz < -halfCarTol) ? q : 0;
405 }
406 else if (xi/(xSemiAxis*xSemiAxis)*v.x()
407 + yi/(ySemiAxis*ySemiAxis)*v.y()
408 {
409 //
410 // Else, if we are traveling outwards,
411 // we must miss
412 //
413 // return kInfinity;
414 }
415 }
416 }
417
418 //
419 // Check z = +dz planer surface
420 //
421 sigz = p.z() - zTopCut;
422
423 if (sigz > -halfCarTol)
424 {
425 if (v.z() >= 0)
426 {
427
428 if (sigz > 0) return kInfinity;
429
430 if ( sqr(p.x()/( xSemiAxis - halfCarTol
431 + sqr(p.y()/( ySemiAxis - halfCarTol
432 return kInfinity;
433
434 }
435 else {
436 G4double q = -sigz/v.z();
437
438 G4double xi = p.x() + q*v.x(),
439 yi = p.y() + q*v.y();
440
441 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxi
442 {
443 return (sigz > -halfCarTol) ? q : 0;
444 }
445 else if (xi/(xSemiAxis*xSemiAxis)*v.x()
446 + yi/(ySemiAxis*ySemiAxis)*v.y()
447 {
448 // return kInfinity;
449 }
450 }
451 }
452
453
454 #if 0
455
456 // check to see if Z plane is relevant
457 //
458 if (p.z() < -zTopCut - halfCarTol)
459 {
460 if (v.z() <= 0.0)
461 return distMin;
462
463 G4double lambda = (-zTopCut - p.z())/v.z()
464
465 if ( sqr((lambda*v.x()+p.x())/xSemiAxis) +
466 sqr((lambda*v.y()+p.y())/ySemiAxis) <
467 sqr(zTopCut + zheight + halfCarTol) )
468 {
469 return distMin = std::fabs(lambda);
470 }
471 }
472
473 if (p.z() > zTopCut + halfCarTol)
474 {
475 if (v.z() >= 0.0)
476 { return distMin; }
477
478 G4double lambda = (zTopCut - p.z()) / v.z
479
480 if ( sqr((lambda*v.x() + p.x())/xSemiAxis)
481 sqr((lambda*v.y() + p.y())/ySemiAxis)
482 sqr(zheight - zTopCut + halfCarTol) )
483 {
484 return distMin = std::fabs(lambda);
485 }
486 }
487
488 if (p.z() > zTopCut - halfCarTol
489 && p.z() < zTopCut + halfCarTol )
490 {
491 if (v.z() > 0.)
492 { return kInfinity; }
493
494 return distMin = 0.;
495 }
496
497 if (p.z() < -zTopCut + halfCarTol
498 && p.z() > -zTopCut - halfCarTol)
499 {
500 if (v.z() < 0.)
501 { return distMin = kInfinity; }
502
503 return distMin = 0.;
504 }
505
506 #endif
507
508 // if we are here then it either intersects
509 // or it does not intersect at all
510 //
511 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y(
512 G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) +
513 v.y()*p.y()/sqr(ySemiAxis) +
514 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y(
515 sqr(zheight - p.z());
516
517 G4double discr = B*B - 4.*A*C;
518
519 // if the discriminant is negative it never
520 //
521 if ( discr < -halfCarTol )
522 { return distMin; }
523
524 // case below is when it hits or grazes the
525 //
526 if ( (discr >= -halfCarTol ) && (discr < hal
527 {
528 return distMin = std::fabs(-B/(2.*A));
529 }
530
531 G4double plus = (-B+std::sqrt(discr))/(2.*A
532 G4double minus = (-B-std::sqrt(discr))/(2.*A
533
534 // Special case::Point on Surface, Check nor
535
536 if ( ( std::fabs(plus) < halfCarTol )||( std
537 {
538 G4ThreeVector truenorm(p.x()/(xSemiAxis*xS
539 p.y()/(ySemiAxis*yS
540 -( p.z() - zheight
541 if ( truenorm*v >= 0) // going outside t
542 {
543 return kInfinity;
544 }
545 else
546 {
547 return 0;
548 }
549 }
550
551 // G4double lambda = std::fabs(plus) < std::
552 G4double lambda = 0;
553
554 if ( minus > halfCarTol && minus < distMin )
555 {
556 lambda = minus ;
557 // check normal vector n * v < 0
558 G4ThreeVector pin = p + lambda*v;
559 if(std::fabs(pin.z())< zTopCut + halfCarTo
560 {
561 G4ThreeVector truenorm(pin.x()/(xSemiAxi
562 pin.y()/(ySemiAxi
563 - ( pin.z() - zhe
564 if ( truenorm*v < 0)
565 { // yes, going inside the solid
566 distMin = lambda;
567 }
568 }
569 }
570 if ( plus > halfCarTol && plus < distMin )
571 {
572 lambda = plus ;
573 // check normal vector n * v < 0
574 G4ThreeVector pin = p + lambda*v;
575 if(std::fabs(pin.z()) < zTopCut + halfCarT
576 {
577 G4ThreeVector truenorm(pin.x()/(xSemiAxi
578 pin.y()/(ySemiAxi
579 - ( pin.z() - zhe
580 if ( truenorm*v < 0)
581 { // yes, going inside the solid
582 distMin = lambda;
583 }
584 }
585 }
586 if (distMin < halfCarTol) distMin=0.;
587 return distMin ;
588 }
589
590 //////////////////////////////////////////////
591 //
592 // Calculate distance (<= actual) to closest s
593 // Return 0 if point inside
594
595 G4double G4EllipticalCone::DistanceToIn(const
596 {
597 G4double hp = std::sqrt(p.x()*p.x()*invXX +
598 G4double ds = (hp - zheight)*cosAxisMin;
599 G4double dz = std::abs(p.z()) - zTopCut;
600 G4double dist = std::max(ds,dz);
601 return (dist > 0) ? dist : 0.;
602 }
603
604 //////////////////////////////////////////////
605 //
606 // Calculate distance to surface of shape from
607 // allowing for tolerance
608
609 G4double G4EllipticalCone::DistanceToOut(const
610 const
611 const
612
613
614 {
615 G4double distMin, lambda;
616 enum surface_e {kPlaneSurf, kCurvedSurf, kNo
617
618 distMin = kInfinity;
619 surface = kNoSurf;
620
621 if (v.z() < 0.0)
622 {
623 lambda = (-p.z() - zTopCut)/v.z();
624
625 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis
626 sqr((p.y() + lambda*v.y())/ySemiAxis
627 sqr(zheight + zTopCut + halfCarTol)
628 {
629 distMin = std::fabs(lambda);
630
631 if (!calcNorm) { return distMin; }
632 }
633 distMin = std::fabs(lambda);
634 surface = kPlaneSurf;
635 }
636
637 if (v.z() > 0.0)
638 {
639 lambda = (zTopCut - p.z()) / v.z();
640
641 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis
642 + sqr((p.y() + lambda*v.y())/ySemiAxis
643 < (sqr(zheight - zTopCut + halfCarTol))
644 {
645 distMin = std::fabs(lambda);
646 if (!calcNorm) { return distMin; }
647 }
648 distMin = std::fabs(lambda);
649 surface = kPlaneSurf;
650 }
651
652 // if we are here then it either intersects
653 // curved surface...
654 //
655 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y(
656 G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis)
657 v.y()*p.y()/sqr(ySemiAxis)
658 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y(
659 - sqr(zheight - p.z());
660
661 G4double discr = B*B - 4.*A*C;
662
663 if ( discr >= - halfCarTol && discr < halfCa
664 {
665 if(!calcNorm) { return distMin = std::fabs
666 }
667
668 else if ( discr > halfCarTol )
669 {
670 G4double plus = (-B+std::sqrt(discr))/(2.
671 G4double minus = (-B-std::sqrt(discr))/(2.
672
673 if ( plus > halfCarTol && minus > halfCarT
674 {
675 // take the shorter distance
676 //
677 lambda = std::fabs(plus) < std::fabs(m
678 }
679 else
680 {
681 // at least one solution is close to zer
682 // so, take small positive solution or z
683 //
684 lambda = plus > -halfCarTol ? plus : 0
685 }
686
687 if ( std::fabs(lambda) < distMin )
688 {
689 if( std::fabs(lambda) > halfCarTol)
690 {
691 distMin = std::fabs(lambda);
692 surface = kCurvedSurf;
693 }
694 else // Point is On the Surface, Check
695 {
696 G4ThreeVector truenorm(p.x()/(xSemiAxi
697 p.y()/(ySemiAxi
698 -( p.z() - zhei
699 if( truenorm.dot(v) > 0 )
700 {
701 distMin = 0.0;
702 surface = kCurvedSurf;
703 }
704 }
705 }
706 }
707
708 // set normal if requested
709 //
710 if (calcNorm)
711 {
712 if (surface == kNoSurf)
713 {
714 *validNorm = false;
715 }
716 else
717 {
718 *validNorm = true;
719 switch (surface)
720 {
721 case kPlaneSurf:
722 {
723 *n = G4ThreeVector(0.,0.,(v.z() > 0.
724 }
725 break;
726
727 case kCurvedSurf:
728 {
729 G4ThreeVector pexit = p + distMin*v;
730 G4ThreeVector truenorm( pexit.x()/(x
731 pexit.y()/(y
732 -( pexit.z()
733 truenorm /= truenorm.mag();
734 *n= truenorm;
735 }
736 break;
737
738 default: // Should never re
739 DumpInfo();
740 std::ostringstream message;
741 G4long oldprc = message.precision(16
742 message << "Undefined side for valid
743 << G4endl
744 << "Position:" << G4endl
745 << " p.x() = " << p.x()/
746 << " p.y() = " << p.y()/
747 << " p.z() = " << p.z()/
748 << "Direction:" << G4endl
749 << " v.x() = " << v.x()
750 << " v.y() = " << v.y()
751 << " v.z() = " << v.z()
752 << "Proposed distance :" <<
753 << " distMin = " << dis
754 message.precision(oldprc);
755 G4Exception("G4EllipticalCone::Dista
756 "GeomSolids1002", JustWa
757 break;
758 }
759 }
760 }
761
762 if (distMin < halfCarTol) { distMin=0; }
763
764 return distMin;
765 }
766
767 //////////////////////////////////////////////
768 //
769 // Calculate distance (<=actual) to closest su
770
771 G4double G4EllipticalCone::DistanceToOut(const
772 {
773 #ifdef G4SPECSDEBUG
774 if( Inside(p) == kOutside )
775 {
776 std::ostringstream message;
777 G4long oldprc = message.precision(16);
778 message << "Point p is outside (!?) of so
779 << "Position:\n"
780 << " p.x() = " << p.x()/mm <<
781 << " p.y() = " << p.y()/mm <<
782 << " p.z() = " << p.z()/mm <<
783 message.precision(oldprc) ;
784 G4Exception("G4Ellipsoid::DistanceToOut(p
785 JustWarning, message);
786 DumpInfo();
787 }
788 #endif
789 G4double hp = std::sqrt(p.x()*p.x()*invXX +
790 G4double ds = (zheight - hp)*cosAxisMin;
791 G4double dz = zTopCut - std::abs(p.z());
792 G4double dist = std::min(ds,dz);
793 return (dist > 0) ? dist : 0.;
794 }
795
796 //////////////////////////////////////////////
797 //
798 // GetEntityType
799
800 G4GeometryType G4EllipticalCone::GetEntityType
801 {
802 return {"G4EllipticalCone"};
803 }
804
805 //////////////////////////////////////////////
806 //
807 // Make a clone of the object
808
809 G4VSolid* G4EllipticalCone::Clone() const
810 {
811 return new G4EllipticalCone(*this);
812 }
813
814 //////////////////////////////////////////////
815 //
816 // Stream object contents to an output stream
817
818 std::ostream& G4EllipticalCone::StreamInfo( st
819 {
820 G4long oldprc = os.precision(16);
821 os << "-------------------------------------
822 << " *** Dump for solid - " << GetName
823 << " =================================
824 << " Solid type: G4EllipticalCone\n"
825 << " Parameters: \n"
826
827 << " semi-axis x: " << xSemiAxis/mm <<
828 << " semi-axis y: " << ySemiAxis/mm <<
829 << " height z: " << zheight/mm << "
830 << " half length in z: " << zTopCut/m
831 << "-------------------------------------
832 os.precision(oldprc);
833
834 return os;
835 }
836
837 //////////////////////////////////////////////
838 //
839 // Return random point on the surface of the s
840
841 G4ThreeVector G4EllipticalCone::GetPointOnSurf
842 {
843 G4double x0 = xSemiAxis*zheight; // x semi a
844 G4double y0 = ySemiAxis*zheight; // y semi a
845 G4double s0 = G4GeomTools::EllipticConeLater
846 G4double kmin = (zTopCut >= zheight ) ? 0. :
847 G4double kmax = (zTopCut >= zheight ) ? 2. :
848
849 // Set areas (base at -Z, side surface, base
850 //
851 G4double szmin = pi*x0*y0*kmax*kmax;
852 G4double szmax = pi*x0*y0*kmin*kmin;
853 G4double sside = s0*(kmax*kmax - kmin*kmin)
854 G4double ssurf[3] = { szmin, sside, szmax };
855 for (auto i=1; i<3; ++i) { ssurf[i] += ssurf
856
857 // Select surface
858 //
859 G4double select = ssurf[2]*G4UniformRand();
860 G4int k = 2;
861 if (select <= ssurf[1]) k = 1;
862 if (select <= ssurf[0]) k = 0;
863
864 // Pick random point on selected surface
865 //
866 G4ThreeVector p;
867 switch(k)
868 {
869 case 0: // base at -Z, uniform distributio
870 {
871 G4double zh = zheight + zTopCut;
872 G4TwoVector rho = G4RandomPointInEllipse
873 p.set(rho.x(),rho.y(),-zTopCut);
874 break;
875 }
876 case 1: // side surface, uniform distribut
877 {
878 G4double zh = G4RandomRadiusInRing(zheig
879 G4double a = x0;
880 G4double b = y0;
881
882 G4double hh = zheight*zheight;
883 G4double aa = a*a;
884 G4double bb = b*b;
885 G4double R = std::max(a,b);
886 G4double mu_max = R*std::sqrt(hh + R*R);
887
888 G4double x,y;
889 for (auto i=0; i<1000; ++i)
890 {
891 G4double phi = CLHEP::twopi*G4UniformRand();
892 x = std::cos(phi);
893 y = std::sin(phi);
894 G4double xx = x*x;
895 G4double yy = y*y;
896 G4double E = hh + aa*xx + bb*yy;
897 G4double F = (aa-bb)*x*y;
898 G4double G = aa*yy + bb*xx;
899 G4double mu = std::sqrt(E*G - F*F);
900 if (mu_max*G4UniformRand() <= mu) brea
901 }
902 p.set(zh*xSemiAxis*x,zh*ySemiAxis*y,zhei
903 break;
904 }
905 case 2: // base at +Z, uniform distributio
906 {
907 G4double zh = zheight - zTopCut;
908 G4TwoVector rho = G4RandomPointInEllipse
909 p.set(rho.x(),rho.y(),zTopCut);
910 break;
911 }
912 }
913 return p;
914 }
915
916 //////////////////////////////////////////////
917 //
918 // Get cubic volume
919
920 G4double G4EllipticalCone::GetCubicVolume()
921 {
922 if (fCubicVolume == 0.0)
923 {
924 G4double x0 = xSemiAxis*zheight; // x semi
925 G4double y0 = ySemiAxis*zheight; // y semi
926 G4double v0 = CLHEP::pi*x0*y0*zheight/3.;
927 G4double kmin = (zTopCut >= zheight ) ? 0.
928 G4double kmax = (zTopCut >= zheight ) ? 2.
929 fCubicVolume = (kmax - kmin)*(kmax*kmax +
930 }
931 return fCubicVolume;
932 }
933
934 //////////////////////////////////////////////
935 //
936 // Get surface area
937
938 G4double G4EllipticalCone::GetSurfaceArea()
939 {
940 if (fSurfaceArea == 0.0)
941 {
942 G4double x0 = xSemiAxis*zheight; // x semi
943 G4double y0 = ySemiAxis*zheight; // y semi
944 G4double s0 = G4GeomTools::EllipticConeLat
945 G4double kmin = (zTopCut >= zheight ) ? 0.
946 G4double kmax = (zTopCut >= zheight ) ? 2.
947 fSurfaceArea = (kmax - kmin)*(kmax + kmin)
948 + CLHEP::pi*x0*y0*(kmin*kmin
949 }
950 return fSurfaceArea;
951 }
952
953 //////////////////////////////////////////////
954 //
955 // Methods for visualisation
956
957 void G4EllipticalCone::DescribeYourselfTo (G4V
958 {
959 scene.AddSolid(*this);
960 }
961
962 G4VisExtent G4EllipticalCone::GetExtent() cons
963 {
964 // Define the sides of the box into which th
965 //
966 G4ThreeVector pmin,pmax;
967 BoundingLimits(pmin,pmax);
968 return { pmin.x(), pmax.x(), pmin.y(), pmax.
969 }
970
971 G4Polyhedron* G4EllipticalCone::CreatePolyhedr
972 {
973 return new G4PolyhedronEllipticalCone(xSemiA
974 }
975
976 G4Polyhedron* G4EllipticalCone::GetPolyhedron
977 {
978 if ( (fpPolyhedron == nullptr)
979 || fRebuildPolyhedron
980 || (fpPolyhedron->GetNumberOfRotationSteps
981 fpPolyhedron->GetNumberOfRotationSteps
982 {
983 G4AutoLock l(&polyhedronMutex);
984 delete fpPolyhedron;
985 fpPolyhedron = CreatePolyhedron();
986 fRebuildPolyhedron = false;
987 l.unlock();
988 }
989 return fpPolyhedron;
990 }
991
992 #endif // !defined(G4GEOM_USE_UELLIPTICALCONE)
993