| 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 for G4Paraboloid class
27 //
28 // Author : Lukas Lindroos (CERN), July 2007
29 // Revised: Tatiana Nikitina (CERN)
30 // -------------------------------------------
31
32 #include "globals.hh"
33
34 #include "G4Paraboloid.hh"
35
36 #if !(defined(G4GEOM_USE_UPARABOLOID) && defin
37
38 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh"
41
42 #include "meshdefs.hh"
43
44 #include "Randomize.hh"
45
46 #include "G4VGraphicsScene.hh"
47 #include "G4VisExtent.hh"
48
49 #include "G4AutoLock.hh"
50
51 namespace
52 {
53 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE
54 }
55
56 using namespace CLHEP;
57
58 //////////////////////////////////////////////
59 //
60 // constructor - check parameters
61 //
62 G4Paraboloid::G4Paraboloid(const G4String& pNa
63 G4double pDz,
64 G4double pR1,
65 G4double pR2)
66 : G4VSolid(pName)
67 {
68 if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.
69 {
70 std::ostringstream message;
71 message << "Invalid dimensions. Negative I
72 << GetName();
73 G4Exception("G4Paraboloid::G4Paraboloid()"
74 FatalErrorInArgument, message,
75 "Z half-length must be larger
76 }
77
78 r1 = pR1;
79 r2 = pR2;
80 dz = pDz;
81
82 // r1^2 = k1 * (-dz) + k2
83 // r2^2 = k1 * ( dz) + k2
84 // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 +
85 // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) =>
86
87 k1 = (r2 * r2 - r1 * r1) / 2 / dz;
88 k2 = (r2 * r2 + r1 * r1) / 2;
89 }
90
91 //////////////////////////////////////////////
92 //
93 // Fake default constructor - sets only member
94 // for usage restri
95 //
96 G4Paraboloid::G4Paraboloid( __void__& a )
97 : G4VSolid(a), dz(0.), r1(0.), r2(0.), k1(0.
98 {
99 }
100
101 //////////////////////////////////////////////
102 //
103 // Destructor
104 //
105 G4Paraboloid::~G4Paraboloid()
106 {
107 delete fpPolyhedron; fpPolyhedron = nullptr;
108 }
109
110 //////////////////////////////////////////////
111 //
112 // Copy constructor
113 //
114 G4Paraboloid::G4Paraboloid(const G4Paraboloid&
115 : G4VSolid(rhs),
116 fSurfaceArea(rhs.fSurfaceArea), fCubicVolu
117 dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs
118 {
119 }
120
121 //////////////////////////////////////////////
122 //
123 // Assignment operator
124 //
125 G4Paraboloid& G4Paraboloid::operator = (const
126 {
127 // Check assignment to self
128 //
129 if (this == &rhs) { return *this; }
130
131 // Copy base class data
132 //
133 G4VSolid::operator=(rhs);
134
135 // Copy data
136 //
137 fSurfaceArea = rhs.fSurfaceArea; fCubicVolu
138 dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 =
139 fRebuildPolyhedron = false;
140 delete fpPolyhedron; fpPolyhedron = nullptr
141
142 return *this;
143 }
144
145 //////////////////////////////////////////////
146 //
147 // Get bounding box
148 //
149 void G4Paraboloid::BoundingLimits(G4ThreeVecto
150 G4ThreeVecto
151 {
152 pMin.set(-r2,-r2,-dz);
153 pMax.set( r2, r2, dz);
154
155 // Check correctness of the bounding box
156 //
157 if (pMin.x() >= pMax.x() || pMin.y() >= pMax
158 {
159 std::ostringstream message;
160 message << "Bad bounding box (min >= max)
161 << GetName() << " !"
162 << "\npMin = " << pMin
163 << "\npMax = " << pMax;
164 G4Exception("G4Paraboloid::BoundingLimits(
165 JustWarning, message);
166 DumpInfo();
167 }
168 }
169
170 //////////////////////////////////////////////
171 //
172 // Calculate extent under transform and specif
173 //
174 G4bool
175 G4Paraboloid::CalculateExtent(const EAxis pAxi
176 const G4VoxelLim
177 const G4AffineTr
178 G4double&
179 {
180 G4ThreeVector bmin, bmax;
181
182 // Get bounding box
183 BoundingLimits(bmin,bmax);
184
185 // Find extent
186 G4BoundingEnvelope bbox(bmin,bmax);
187 return bbox.CalculateExtent(pAxis,pVoxelLimi
188 }
189
190 //////////////////////////////////////////////
191 //
192 // Return whether point inside/outside/on surf
193 //
194 EInside G4Paraboloid::Inside(const G4ThreeVect
195 {
196 // First check is the point is above or bel
197 //
198 if(std::fabs(p.z()) > dz + 0.5 * kCarToleran
199
200 G4double rho2 = p.perp2(),
201 rhoSurfTimesTol2 = (k1 * p.z() + k
202 A = rho2 - ((k1 *p.z() + k2) + 0.25
203
204 if(A < 0 && sqr(A) > rhoSurfTimesTol2)
205 {
206 // Actually checking rho < radius of parab
207 // We're either inside or in lower/upper c
208
209 if(std::fabs(p.z()) > dz - 0.5 * kCarToler
210 {
211 // We're in the upper/lower cutoff area,
212 // maybe further checks should be made t
213
214 return kSurface;
215 }
216 else
217 {
218 return kInside;
219 }
220 }
221 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
222 {
223 // We're in the parabolic surface.
224
225 return kSurface;
226 }
227 else
228 {
229 return kOutside;
230 }
231 }
232
233 //////////////////////////////////////////////
234 //
235 // SurfaceNormal
236 //
237 G4ThreeVector G4Paraboloid::SurfaceNormal( con
238 {
239 G4ThreeVector n(0, 0, 0);
240 if(std::fabs(p.z()) > dz + 0.5 * kCarToleran
241 {
242 // If above or below just return normal ve
243
244 n = G4ThreeVector(0, 0, p.z()/std::fabs(p.
245 }
246 else if(std::fabs(p.z()) > dz - 0.5 * kCarTo
247 {
248 // This means we're somewhere in the plane
249 // (As far as the program is concerned any
250
251 if(p.z() < 0) // Are we in upper or lower
252 {
253 if(p.perp2() > sqr(r1 + 0.5 * kCarTolera
254 {
255 n = G4ThreeVector(p.x(), p.y(), -k1 /
256 }
257 else if(r1 < 0.5 * kCarTolerance
258 || p.perp2() > sqr(r1 - 0.5 * kCarT
259 {
260 n = G4ThreeVector(p.x(), p.y(), 0.).un
261 + G4ThreeVector(0., 0., -1.).unit();
262 n = n.unit();
263 }
264 else
265 {
266 n = G4ThreeVector(0., 0., -1.);
267 }
268 }
269 else
270 {
271 if(p.perp2() > sqr(r2 + 0.5 * kCarTolera
272 {
273 n = G4ThreeVector(p.x(), p.y(), 0.).un
274 }
275 else if(r2 < 0.5 * kCarTolerance
276 || p.perp2() > sqr(r2 - 0.5 * kCarT
277 {
278 n = G4ThreeVector(p.x(), p.y(), 0.).un
279 + G4ThreeVector(0., 0., 1.).unit();
280 n = n.unit();
281 }
282 else
283 {
284 n = G4ThreeVector(0., 0., 1.);
285 }
286 }
287 }
288 else
289 {
290 G4double rho2 = p.perp2();
291 G4double rhoSurfTimesTol2 = (k1 * p.z() +
292 G4double A = rho2 - ((k1 *p.z() + k2)
293 + 0.25 * kCarTolerance * kCarTo
294
295 if(A < 0 && sqr(A) > rhoSurfTimesTol2)
296 {
297 // Actually checking rho < radius of par
298 // We're inside.
299
300 if(p.mag2() != 0) { n = p.unit(); }
301 }
302 else if(A <= 0 || sqr(A) < rhoSurfTimesTol
303 {
304 // We're in the parabolic surface.
305
306 n = G4ThreeVector(p.x(), p.y(), - k1 / 2
307 }
308 else
309 {
310 n = G4ThreeVector(p.x(), p.y(), - k1 / 2
311 }
312 }
313
314 if(n.mag2() == 0)
315 {
316 std::ostringstream message;
317 message << "No normal defined for this poi
318 << " p = " << 1 / mm * p
319 G4Exception("G4Paraboloid::SurfaceNormal(p
320 JustWarning, message);
321 }
322 return n;
323 }
324
325 //////////////////////////////////////////////
326 //
327 // Calculate distance to shape from outside, a
328 // - return kInfinity if no intersection
329 //
330 //
331 G4double G4Paraboloid::DistanceToIn( const G4T
332 const G4T
333 {
334 G4double rho2 = p.perp2(), paraRho2 = std::f
335 G4double tol2 = kCarTolerance*kCarTolerance;
336 G4double tolh = 0.5*kCarTolerance;
337
338 if((r2 != 0.0) && p.z() > - tolh + dz)
339 {
340 // If the points is above check for inters
341
342 if(v.z() < 0)
343 {
344 G4double intersection = (dz - p.z()) / v
345 if(sqr(p.x() + v.x()*intersection)
346 + sqr(p.y() + v.y()*intersection) < sqr
347 {
348 if(p.z() < tolh + dz)
349 { return 0; }
350 else
351 { return intersection; }
352 }
353 }
354 else // Direction away, no possibility of
355 {
356 return kInfinity;
357 }
358 }
359 else if((r1 != 0.0) && p.z() < tolh - dz)
360 {
361 // If the points is belove check for inter
362
363 if(v.z() > 0)
364 {
365 G4double intersection = (-dz - p.z()) /
366 if(sqr(p.x() + v.x()*intersection)
367 + sqr(p.y() + v.y()*intersection) < sqr
368 {
369 if(p.z() > -tolh - dz)
370 {
371 return 0;
372 }
373 else
374 {
375 return intersection;
376 }
377 }
378 }
379 else // Direction away, no possibility of
380 {
381 return kInfinity;
382 }
383 }
384
385 G4double A = k1 / 2 * v.z() - p.x() * v.x()
386 vRho2 = v.perp2(), intersection,
387 B = (k1 * p.z() + k2 - rho2) * vRho
388
389 if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRh
390 || (p.z() < - dz+kCarTolerance)
391 || (p.z() > dz-kCarTolerance) ) // Make su
392 {
393 // Is there a problem with squaring rho tw
394
395 if(vRho2<tol2) // Needs to be treated sepe
396 {
397 intersection = ((rho2 - k2)/k1 - p.z())/
398 if(intersection < 0) { return kInfinity;
399 else if(std::fabs(p.z() + v.z() * inters
400 {
401 return intersection;
402 }
403 else
404 {
405 return kInfinity;
406 }
407 }
408 else if(A*A + B < 0) // No real intersecti
409 {
410 return kInfinity;
411 }
412 else
413 {
414 intersection = (A - std::sqrt(B + sqr(A)
415 if(intersection < 0)
416 {
417 return kInfinity;
418 }
419 else if(std::fabs(p.z() + intersection *
420 {
421 return intersection;
422 }
423 else
424 {
425 return kInfinity;
426 }
427 }
428 }
429 else if(sqr(rho2 - paraRho2 - .25 * tol2) <=
430 {
431 // If this is true we're somewhere in the
432
433 G4ThreeVector normal(p.x(), p.y(), -k1/2);
434 if(normal.dot(v) <= 0)
435 { return 0; }
436 else
437 { return kInfinity; }
438 }
439 else
440 {
441 std::ostringstream message;
442 if(Inside(p) == kInside)
443 {
444 message << "Point p is inside! - " << Ge
445 }
446 else
447 {
448 message << "Likely a problem in this fun
449 << G4endl;
450 }
451 message << " p = " << p * (1/mm)
452 << " v = " << v * (1/mm)
453 G4Exception("G4Paraboloid::DistanceToIn(p,
454 JustWarning, message);
455 return 0;
456 }
457 }
458
459 //////////////////////////////////////////////
460 //
461 // Calculate distance (<= actual) to closest s
462 // - Return zero if point inside
463 //
464 G4double G4Paraboloid::DistanceToIn(const G4Th
465 {
466 G4double safz = -dz+std::fabs(p.z());
467 if(safz<0.) { safz=0.; }
468 G4double safr = kInfinity;
469
470 G4double rho = p.x()*p.x()+p.y()*p.y();
471 G4double paraRho = (p.z()-k2)/k1;
472 G4double sqrho = std::sqrt(rho);
473
474 if(paraRho<0.)
475 {
476 safr=sqrho-r2;
477 if(safr>safz) { safz=safr; }
478 return safz;
479 }
480
481 G4double sqprho = std::sqrt(paraRho);
482 G4double dRho = sqrho-sqprho;
483 if(dRho<0.) { return safz; }
484
485 G4double talf = -2.*k1*sqprho;
486 G4double tmp = 1+talf*talf;
487 if(tmp<0.) { return safz; }
488
489 G4double salf = talf/std::sqrt(tmp);
490 safr = std::fabs(dRho*salf);
491 if(safr>safz) { safz=safr; }
492
493 return safz;
494 }
495
496 //////////////////////////////////////////////
497 //
498 // Calculate distance to surface of shape from
499 //
500 G4double G4Paraboloid::DistanceToOut(const G4T
501 const G4Th
502 const G4bo
503 G4bo
504 G4Th
505 {
506 G4double rho2 = p.perp2(), paraRho2 = std::f
507 G4double vRho2 = v.perp2(), intersection;
508 G4double tol2 = kCarTolerance*kCarTolerance;
509 G4double tolh = 0.5*kCarTolerance;
510
511 if(calcNorm) { *validNorm = false; }
512
513 // We have that the particle p follows the l
514 // meaning x = p.x() + s * v.x(), y = p.y()
515 // z = p.z() + s * v.z()
516 // The equation for all points on the surfac
517 // to include all z) x^2 + y^2 = k1 * z + k2
518 // => s = (A +- std::sqrt(A^2 + B)) / vRho2
519 // where:
520 //
521 G4double A = k1 / 2 * v.z() - p.x() * v.x()
522 //
523 // and:
524 //
525 G4double B = (-rho2 + paraRho2) * vRho2;
526
527 if ( rho2 < paraRho2 && sqr(rho2 - paraRho2
528 && std::fabs(p.z()) < dz - kCarTolerance)
529 {
530 // Make sure it's safely inside.
531
532 if(v.z() > 0)
533 {
534 // It's heading upwards, check where it
535 // When it does, is that in the surface
536 // z = p.z() + variable * v.z() for all
537 // => variable = (z - p.z()) / v.z() so
538
539 intersection = (dz - p.z()) / v.z();
540 G4ThreeVector ip = p + intersection * v;
541
542 if(ip.perp2() < sqr(r2 + kCarTolerance))
543 {
544 if(calcNorm)
545 {
546 *n = G4ThreeVector(0, 0, 1);
547 if(r2 < tolh || ip.perp2() > sqr(r2
548 {
549 *n += G4ThreeVector(ip.x(), ip.y()
550 *n = n->unit();
551 }
552 *validNorm = true;
553 }
554 return intersection;
555 }
556 }
557 else if(v.z() < 0)
558 {
559 // It's heading downwards, check were it
560 // When it does, is that in the surface
561 // z = p.z() + variable * v.z() for all
562 // => variable = (z - p.z()) / v.z() so
563
564 intersection = (-dz - p.z()) / v.z();
565 G4ThreeVector ip = p + intersection * v;
566
567 if(ip.perp2() < sqr(r1 + tolh))
568 {
569 if(calcNorm)
570 {
571 *n = G4ThreeVector(0, 0, -1);
572 if(r1 < tolh || ip.perp2() > sqr(r1
573 {
574 *n += G4ThreeVector(ip.x(), ip.y()
575 *n = n->unit();
576 }
577 *validNorm = true;
578 }
579 return intersection;
580 }
581 }
582
583 // Now check for collisions with paraboloi
584
585 if(vRho2 == 0) // Needs to be treated sepe
586 {
587 intersection = ((rho2 - k2)/k1 - p.z())/
588 if(calcNorm)
589 {
590 G4ThreeVector intersectionP = p + v *
591 *n = G4ThreeVector(intersectionP.x(),
592 *n = n->unit();
593
594 *validNorm = true;
595 }
596 return intersection;
597 }
598 else if( ((A <= 0) && (B >= sqr(A) * (sqr(
599 {
600 // intersection = (A + std::sqrt(B + sqr
601 // The above calculation has a precision
602 // known problem of solving quadratic eq
603
604 A = A/vRho2;
605 B = (k1 * p.z() + k2 - rho2)/vRho2;
606 intersection = B/(-A + std::sqrt(B + sqr
607 if(calcNorm)
608 {
609 G4ThreeVector intersectionP = p + v *
610 *n = G4ThreeVector(intersectionP.x(),
611 *n = n->unit();
612 *validNorm = true;
613 }
614 return intersection;
615 }
616 std::ostringstream message;
617 message << "There is no intersection betwe
618 << G4endl
619 << " p = " << p << G4endl
620 << " v = " << v;
621 G4Exception("G4Paraboloid::DistanceToOut(p
622 JustWarning, message);
623
624 return kInfinity;
625 }
626 else if ( (rho2 < paraRho2 + kCarTolerance
627 || sqr(rho2 - paraRho2 - 0.25 * tol2)
628 && std::fabs(p.z()) < dz + tolh)
629 {
630 // If this is true we're somewhere in the
631
632 G4ThreeVector normal = G4ThreeVector (p.x(
633
634 if(std::fabs(p.z()) > dz - tolh)
635 {
636 // We're in the lower or upper edge
637 //
638 if( ((v.z() > 0) && (p.z() > 0)) || ((v.
639 { // If we're heading out of
640 if(calcNorm)
641 {
642 *validNorm = true;
643 if(p.z() > 0)
644 { *n = G4ThreeVector(0, 0, 1); }
645 else
646 { *n = G4ThreeVector(0, 0, -1); }
647 }
648 return 0;
649 }
650
651 if(v.z() == 0)
652 {
653 // Case where we're moving inside the
654 // treated separately.
655 // Distance until it goes out through
656
657 G4double r = (p.z() > 0)? r2 : r1;
658 G4double pDotV = p.dot(v);
659 A = vRho2 * ( sqr(r) - sqr(p.x()) - sq
660 intersection = (-pDotV + std::sqrt(A +
661
662 if(calcNorm)
663 {
664 *validNorm = true;
665
666 *n = (G4ThreeVector(0, 0, p.z()/std:
667 + G4ThreeVector(p.x() + v.x() *
668 * intersection, -k1/2).unit()).u
669 }
670 return intersection;
671 }
672 }
673 //
674 // Problem in the Logic :: Following condi
675 // and Vz<0 will
676 // it has to retur
677 // surface or with
678 // The logic has to be :: If not found in
679 // do not exit but continue to search for
680 // Only for point situated on both borders
681 // this condition has to be taken into acc
682 //
683 //
684 // else if(normal.dot(v) >= 0)
685 // {
686 // if(calcNorm)
687 // {
688 // *validNorm = true;
689 // *n = normal.unit();
690 // }
691 // return 0;
692 // }
693
694 if(v.z() > 0)
695 {
696 // Check for collision with upper edge.
697
698 intersection = (dz - p.z()) / v.z();
699 G4ThreeVector ip = p + intersection * v;
700
701 if(ip.perp2() < sqr(r2 - tolh))
702 {
703 if(calcNorm)
704 {
705 *validNorm = true;
706 *n = G4ThreeVector(0, 0, 1);
707 }
708 return intersection;
709 }
710 else if(ip.perp2() < sqr(r2 + tolh))
711 {
712 if(calcNorm)
713 {
714 *validNorm = true;
715 *n = G4ThreeVector(0, 0, 1)
716 + G4ThreeVector(ip.x(), ip.y(), -
717 *n = n->unit();
718 }
719 return intersection;
720 }
721 }
722 if( v.z() < 0)
723 {
724 // Check for collision with lower edge.
725
726 intersection = (-dz - p.z()) / v.z();
727 G4ThreeVector ip = p + intersection * v;
728
729 if(ip.perp2() < sqr(r1 - tolh))
730 {
731 if(calcNorm)
732 {
733 *validNorm = true;
734 *n = G4ThreeVector(0, 0, -1);
735 }
736 return intersection;
737 }
738 else if(ip.perp2() < sqr(r1 + tolh))
739 {
740 if(calcNorm)
741 {
742 *validNorm = true;
743 *n = G4ThreeVector(0, 0, -1)
744 + G4ThreeVector(ip.x(), ip.y(), -
745 *n = n->unit();
746 }
747 return intersection;
748 }
749 }
750
751 // Note: comparison with zero below would
752 //
753 if(std::fabs(vRho2) > tol2) // precision e
754 { // intersectio
755 A = A/vRho2;
756 B = (k1 * p.z() + k2 - rho2);
757 if(std::fabs(B)>kCarTolerance)
758 {
759 B = (B)/vRho2;
760 intersection = B/(-A + std::sqrt(B + s
761 }
762 else // Point is On
763 { // solution de
764 if(normal.dot(v) >= 0)
765 {
766 if(calcNorm)
767 {
768 *validNorm = true;
769 *n = normal.unit();
770 }
771 return 0;
772 }
773 intersection = 2.*A;
774 }
775 }
776 else
777 {
778 intersection = ((rho2 - k2) / k1 - p.z()
779 }
780
781 if(calcNorm)
782 {
783 *validNorm = true;
784 *n = G4ThreeVector(p.x() + intersection
785 + intersection * v.y(), - k1 / 2);
786 *n = n->unit();
787 }
788 return intersection;
789 }
790 else
791 {
792 #ifdef G4SPECSDEBUG
793 if(kOutside == Inside(p))
794 {
795 G4Exception("G4Paraboloid::DistanceToOut
796 JustWarning, "Point p is out
797 }
798 else
799 G4Exception("G4Paraboloid::DistanceToOut
800 JustWarning, "There's an err
801 #endif
802 return kInfinity;
803 }
804 return 0;
805 }
806
807 //////////////////////////////////////////////
808 //
809 // Calculate distance (<=actual) to closest su
810 //
811 G4double G4Paraboloid::DistanceToOut(const G4T
812 {
813 G4double safe=0.0,rho,safeR,safeZ ;
814 G4double tanRMax,secRMax,pRMax ;
815
816 #ifdef G4SPECSDEBUG
817 if( Inside(p) == kOutside )
818 {
819 G4cout << G4endl ;
820 DumpInfo();
821 std::ostringstream message;
822 G4long oldprc = message.precision(16);
823 message << "Point p is outside !?" << G4e
824 << "Position:" << G4endl
825 << " p.x() = " << p.x()/mm <<
826 << " p.y() = " << p.y()/mm <<
827 << " p.z() = " << p.z()/mm <<
828 message.precision(oldprc) ;
829 G4Exception("G4Paraboloid::DistanceToOut(
830 JustWarning, message);
831 }
832 #endif
833
834 rho = p.perp();
835 safeZ = dz - std::fabs(p.z()) ;
836
837 tanRMax = (r2 - r1)*0.5/dz ;
838 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
839 pRMax = tanRMax*p.z() + (r1+r2)*0.5 ;
840 safeR = (pRMax - rho)/secRMax ;
841
842 if (safeZ < safeR) { safe = safeZ; }
843 else { safe = safeR; }
844 if ( safe < 0.5 * kCarTolerance ) { safe = 0
845 return safe ;
846 }
847
848 //////////////////////////////////////////////
849 //
850 // G4EntityType
851 //
852 G4GeometryType G4Paraboloid::GetEntityType() c
853 {
854 return {"G4Paraboloid"};
855 }
856
857 //////////////////////////////////////////////
858 //
859 // Make a clone of the object
860 //
861 G4VSolid* G4Paraboloid::Clone() const
862 {
863 return new G4Paraboloid(*this);
864 }
865
866 //////////////////////////////////////////////
867 //
868 // Stream object contents to an output stream
869 //
870 std::ostream& G4Paraboloid::StreamInfo( std::o
871 {
872 G4long oldprc = os.precision(16);
873 os << "-------------------------------------
874 << " *** Dump for solid - " << GetName
875 << " =================================
876 << " Solid type: G4Paraboloid\n"
877 << " Parameters: \n"
878 << " z half-axis: " << dz/mm << " mm
879 << " radius at -dz: " << r1/mm << " mm
880 << " radius at dz: " << r2/mm << " mm
881 << "-------------------------------------
882 os.precision(oldprc);
883
884 return os;
885 }
886
887 //////////////////////////////////////////////
888 //
889 // GetPointOnSurface
890 //
891 G4ThreeVector G4Paraboloid::GetPointOnSurface(
892 {
893 G4double A = (fSurfaceArea == 0)? CalculateS
894 G4double z = G4RandFlat::shoot(0.,1.);
895 G4double phi = G4RandFlat::shoot(0., twopi);
896 if(pi*(sqr(r1) + sqr(r2))/A >= z)
897 {
898 G4double rho;
899 if(pi * sqr(r1) / A > z)
900 {
901 rho = r1 * std::sqrt(G4RandFlat::shoot(0
902 return { rho * std::cos(phi), rho * std:
903 }
904 else
905 {
906 rho = r2 * std::sqrt(G4RandFlat::shoot(0
907 return { rho * std::cos(phi), rho * std:
908 }
909 }
910 else
911 {
912 z = G4RandFlat::shoot(0., 1.)*2*dz - dz;
913 return { std::sqrt(z*k1 + k2)*std::cos(phi
914 std::sqrt(z*k1 + k2)*std::sin(phi
915 }
916 }
917
918 //////////////////////////////////////////////
919 //
920 // Methods for visualisation
921 //
922 void G4Paraboloid::DescribeYourselfTo (G4VGrap
923 {
924 scene.AddSolid(*this);
925 }
926
927 G4Polyhedron* G4Paraboloid::CreatePolyhedron (
928 {
929 return new G4PolyhedronParaboloid(r1, r2, dz
930 }
931
932 G4Polyhedron* G4Paraboloid::GetPolyhedron () c
933 {
934 if (fpPolyhedron == nullptr ||
935 fRebuildPolyhedron ||
936 fpPolyhedron->GetNumberOfRotationStepsAt
937 fpPolyhedron->GetNumberOfRotationSteps()
938 {
939 G4AutoLock l(&polyhedronMutex);
940 delete fpPolyhedron;
941 fpPolyhedron = CreatePolyhedron();
942 fRebuildPolyhedron = false;
943 l.unlock();
944 }
945 return fpPolyhedron;
946 }
947
948 #endif // !defined(G4GEOM_USE_UPARABOLOID) ||
949