| 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 // class G4Ellipsoid
27 //
28 // Implementation of G4Ellipsoid class
29 //
30 // 10.11.99 G.Horton-Smith: first writing, bas
31 // 25.02.05 G.Guerrieri: Revised
32 // 15.12.19 E.Tcherniaev: Complete revision
33 // -------------------------------------------
34
35 #include "G4Ellipsoid.hh"
36
37 #if !(defined(G4GEOM_USE_UELLIPSOID) && define
38
39 #include "globals.hh"
40
41 #include "G4VoxelLimits.hh"
42 #include "G4AffineTransform.hh"
43 #include "G4GeometryTolerance.hh"
44 #include "G4BoundingEnvelope.hh"
45 #include "G4RandomTools.hh"
46 #include "G4QuickRand.hh"
47
48 #include "G4VPVParameterisation.hh"
49
50 #include "G4VGraphicsScene.hh"
51 #include "G4VisExtent.hh"
52
53 #include "G4AutoLock.hh"
54
55 namespace
56 {
57 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZ
58 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZ
59 }
60
61 using namespace CLHEP;
62
63 //////////////////////////////////////////////
64 //
65 // Constructor
66
67 G4Ellipsoid::G4Ellipsoid(const G4String& name,
68 G4double xSemiA
69 G4double ySemiA
70 G4double zSemiA
71 G4double zBotto
72 G4double zTopCu
73 : G4VSolid(name), fDx(xSemiAxis), fDy(ySemiA
74 fZBottomCut(zBottomCut), fZTopCut(zTopCut)
75 {
76 CheckParameters();
77 }
78
79 //////////////////////////////////////////////
80 //
81 // Fake default constructor - sets only member
82 // for usage restri
83
84 G4Ellipsoid::G4Ellipsoid( __void__& a )
85 : G4VSolid(a), fDx(0.), fDy(0.), fDz(0.), fZ
86 {
87 }
88
89 //////////////////////////////////////////////
90 //
91 // Destructor
92
93 G4Ellipsoid::~G4Ellipsoid()
94 {
95 delete fpPolyhedron; fpPolyhedron = nullptr;
96 }
97
98 //////////////////////////////////////////////
99 //
100 // Copy constructor
101
102 G4Ellipsoid::G4Ellipsoid(const G4Ellipsoid& rh
103 : G4VSolid(rhs),
104 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
105 fZBottomCut(rhs.fZBottomCut), fZTopCut(rhs.
106 halfTolerance(rhs.halfTolerance),
107 fXmax(rhs.fXmax), fYmax(rhs.fYmax),
108 fRsph(rhs.fRsph), fR(rhs.fR),
109 fSx(rhs.fSx), fSy(rhs.fSy), fSz(rhs.fSz),
110 fZMidCut(rhs.fZMidCut), fZDimCut(rhs.fZDimC
111 fQ1(rhs.fQ1), fQ2(rhs.fQ2),
112 fCubicVolume(rhs.fCubicVolume),
113 fSurfaceArea(rhs.fSurfaceArea),
114 fLateralArea(rhs.fLateralArea)
115 {
116 }
117
118 //////////////////////////////////////////////
119 //
120 // Assignment operator
121
122 G4Ellipsoid& G4Ellipsoid::operator = (const G4
123 {
124 // Check assignment to self
125 //
126 if (this == &rhs) { return *this; }
127
128 // Copy base class data
129 //
130 G4VSolid::operator=(rhs);
131
132 // Copy data
133 //
134 fDx = rhs.fDx;
135 fDy = rhs.fDy;
136 fDz = rhs.fDz;
137 fZBottomCut = rhs.fZBottomCut;
138 fZTopCut = rhs.fZTopCut;
139
140 halfTolerance = rhs.halfTolerance;
141 fXmax = rhs.fXmax;
142 fYmax = rhs.fYmax;
143 fRsph = rhs.fRsph;
144 fR = rhs.fR;
145 fSx = rhs.fSx;
146 fSy = rhs.fSy;
147 fSz = rhs.fSz;
148 fZMidCut = rhs.fZMidCut;
149 fZDimCut = rhs.fZDimCut;
150 fQ1 = rhs.fQ1;
151 fQ2 = rhs.fQ2;
152
153 fCubicVolume = rhs.fCubicVolume;
154 fSurfaceArea = rhs.fSurfaceArea;
155 fLateralArea = rhs.fLateralArea;
156
157 fRebuildPolyhedron = false;
158 delete fpPolyhedron; fpPolyhedron = nullptr
159
160 return *this;
161 }
162
163 //////////////////////////////////////////////
164 //
165 // Check parameters and make precalculation
166
167 void G4Ellipsoid::CheckParameters()
168 {
169 halfTolerance = 0.5 * kCarTolerance; // half
170 G4double dmin = 2. * kCarTolerance;
171
172 // Check dimensions
173 //
174 if (fDx < dmin || fDy < dmin || fDz < dmin)
175 {
176 std::ostringstream message;
177 message << "Invalid (too small or negative
178 << GetName() << "\n"
179 << " semi-axis x: " << fDx << "\n
180 << " semi-axis y: " << fDy << "\n
181 << " semi-axis z: " << fDz;
182 G4Exception("G4Ellipsoid::CheckParameters(
183 FatalException, message);
184 }
185 G4double A = fDx;
186 G4double B = fDy;
187 G4double C = fDz;
188
189 // Check cuts
190 //
191 if (fZBottomCut == 0. && fZTopCut == 0.)
192 {
193 fZBottomCut = -C;
194 fZTopCut = C;
195 }
196 if (fZBottomCut >= C || fZTopCut <= -C || fZ
197 {
198 std::ostringstream message;
199 message << "Invalid Z cuts for Solid: "
200 << GetName() << "\n"
201 << " bottom cut: " << fZBottomCut
202 << " top cut: " << fZTopCut;
203 G4Exception("G4Ellipsoid::CheckParameters(
204 FatalException, message);
205
206 }
207 fZBottomCut = std::max(fZBottomCut, -C);
208 fZTopCut = std::min(fZTopCut, C);
209
210 // Set extent in x and y
211 fXmax = A;
212 fYmax = B;
213 if (fZBottomCut > 0.)
214 {
215 G4double ratio = fZBottomCut / C;
216 G4double scale = std::sqrt((1. - ratio) *
217 fXmax *= scale;
218 fYmax *= scale;
219 }
220 if (fZTopCut < 0.)
221 {
222 G4double ratio = fZTopCut / C;
223 G4double scale = std::sqrt((1. - ratio) *
224 fXmax *= scale;
225 fYmax *= scale;
226 }
227
228 // Set scale factors
229 fRsph = std::max(std::max(A, B), C); // boun
230 fR = std::min(std::min(A, B), C); // radi
231 fSx = fR / A; // X scale factor
232 fSy = fR / B; // Y scale factor
233 fSz = fR / C; // Z scale factor
234
235 // Scaled cuts
236 fZMidCut = 0.5 * (fZTopCut + fZBottomCut) *
237 fZDimCut = 0.5 * (fZTopCut - fZBottomCut) *
238
239 // Coefficients for approximation of distanc
240 fQ1 = 0.5 / fR;
241 fQ2 = 0.5 * fR + halfTolerance * halfToleran
242
243 fCubicVolume = 0.; // volume
244 fSurfaceArea = 0.; // surface area
245 fLateralArea = 0.; // lateral surface area
246 }
247
248 //////////////////////////////////////////////
249 //
250 // Dispatch to parameterisation for replicatio
251 // computation & modification
252
253 void G4Ellipsoid::ComputeDimensions(G4VPVParam
254 const G4in
255 const G4VP
256 {
257 p->ComputeDimensions(*this,n,pRep);
258 }
259
260 //////////////////////////////////////////////
261 //
262 // Get bounding box
263
264 void G4Ellipsoid::BoundingLimits(G4ThreeVector
265 G4ThreeVector
266 {
267 pMin.set(-fXmax,-fYmax, fZBottomCut);
268 pMax.set( fXmax, fYmax, fZTopCut);
269 }
270
271 //////////////////////////////////////////////
272 //
273 // Calculate extent under transform and specif
274
275 G4bool
276 G4Ellipsoid::CalculateExtent(const EAxis pAxis
277 const G4VoxelLimi
278 const G4AffineTra
279 G4double& p
280 {
281 G4ThreeVector bmin, bmax;
282
283 // Get bounding box
284 BoundingLimits(bmin,bmax);
285
286 // Find extent
287 G4BoundingEnvelope bbox(bmin,bmax);
288 return bbox.CalculateExtent(pAxis,pVoxelLimi
289 }
290
291 //////////////////////////////////////////////
292 //
293 // Return position of point: inside/outside/on
294
295 EInside G4Ellipsoid::Inside(const G4ThreeVecto
296 {
297 G4double x = p.x() * fSx;
298 G4double y = p.y() * fSy;
299 G4double z = p.z() * fSz;
300 G4double rr = x * x + y * y + z * z;
301 G4double distZ = std::abs(z - fZMidCut) - fZ
302 G4double distR = fQ1 * rr - fQ2;
303 G4double dist = std::max(distZ, distR);
304
305 if (dist > halfTolerance) return kOutside;
306 return (dist > -halfTolerance) ? kSurface :
307 }
308
309 //////////////////////////////////////////////
310 //
311 // Return unit normal to surface at p
312
313 G4ThreeVector G4Ellipsoid::SurfaceNormal( cons
314 {
315 G4ThreeVector norm(0., 0., 0.);
316 G4int nsurf = 0;
317
318 // Check cuts
319 G4double x = p.x() * fSx;
320 G4double y = p.y() * fSy;
321 G4double z = p.z() * fSz;
322 G4double distZ = std::abs(z - fZMidCut) - fZ
323 if (std::abs(distZ) <= halfTolerance)
324 {
325 norm.setZ(std::copysign(1., z - fZMidCut))
326 ++nsurf;
327 }
328
329 // Check lateral surface
330 G4double distR = fQ1*(x*x + y*y + z*z) - fQ2
331 if (std::abs(distR) <= halfTolerance)
332 {
333 // normal = (p.x/a^2, p.y/b^2, p.z/c^2)
334 norm += G4ThreeVector(x*fSx, y*fSy, z*fSz)
335 ++nsurf;
336 }
337
338 // Return normal
339 if (nsurf == 1) return norm;
340 else if (nsurf > 1) return norm.unit(); // e
341 else
342 {
343 #ifdef G4SPECSDEBUG
344 std::ostringstream message;
345 G4long oldprc = message.precision(16);
346 message << "Point p is not on surface (!?)
347 << GetName() << "\n";
348 message << "Position:\n";
349 message << " p.x() = " << p.x()/mm << "
350 message << " p.y() = " << p.y()/mm << "
351 message << " p.z() = " << p.z()/mm << "
352 G4cout.precision(oldprc);
353 G4Exception("G4Ellipsoid::SurfaceNormal(p)
354 JustWarning, message );
355 DumpInfo();
356 #endif
357 return ApproxSurfaceNormal(p);
358 }
359 }
360
361 //////////////////////////////////////////////
362 //
363 // Find surface nearest to point and return co
364 // This method normally should not be called.
365
366 G4ThreeVector G4Ellipsoid::ApproxSurfaceNormal
367 {
368 G4double x = p.x() * fSx;
369 G4double y = p.y() * fSy;
370 G4double z = p.z() * fSz;
371 G4double rr = x * x + y * y + z * z;
372 G4double distZ = std::abs(z - fZMidCut) - fZ
373 G4double distR = std::sqrt(rr) - fR;
374 if (distR > distZ && rr > 0.) // distR > di
375 return G4ThreeVector(x*fSx, y*fSy, z*fSz).
376 else
377 return { 0., 0., std::copysign(1., z - fZM
378 }
379
380 //////////////////////////////////////////////
381 //
382 // Calculate distance to shape from outside al
383
384 G4double G4Ellipsoid::DistanceToIn(const G4Thr
385 const G4Thr
386 {
387 G4double offset = 0.;
388 G4ThreeVector pcur = p;
389
390 // Check if point is flying away, relative t
391 //
392 G4double safex = std::abs(p.x()) - fXmax;
393 G4double safey = std::abs(p.y()) - fYmax;
394 G4double safet = p.z() - fZTopCut;
395 G4double safeb = fZBottomCut - p.z();
396
397 if (safex >= -halfTolerance && p.x() * v.x()
398 if (safey >= -halfTolerance && p.y() * v.y()
399 if (safet >= -halfTolerance && v.z() >= 0.)
400 if (safeb >= -halfTolerance && v.z() <= 0.)
401
402 // Relocate point, if required
403 //
404 G4double safe = std::max(std::max(std::max(s
405 if (safe > 32. * fRsph)
406 {
407 offset = (1. - 1.e-08) * safe - 2. * fRsph
408 pcur += offset * v;
409 G4double dist = DistanceToIn(pcur, v);
410 return (dist == kInfinity) ? kInfinity : d
411 }
412
413 // Scale ellipsoid to sphere
414 //
415 G4double px = pcur.x() * fSx;
416 G4double py = pcur.y() * fSy;
417 G4double pz = pcur.z() * fSz;
418 G4double vx = v.x() * fSx;
419 G4double vy = v.y() * fSy;
420 G4double vz = v.z() * fSz;
421
422 // Check if point is leaving the solid
423 //
424 G4double dzcut = fZDimCut;
425 G4double pzcut = pz - fZMidCut;
426 G4double distZ = std::abs(pzcut) - dzcut;
427 if (distZ >= -halfTolerance && pzcut * vz >=
428
429 G4double rr = px * px + py * py + pz * pz;
430 G4double pv = px * vx + py * vy + pz * vz;
431 G4double distR = fQ1 * rr - fQ2;
432 if (distR >= -halfTolerance && pv >= 0.) ret
433
434 G4double A = vx * vx + vy * vy + vz * vz;
435 G4double B = pv;
436 G4double C = rr - fR * fR;
437 G4double D = B * B - A * C;
438 // scratch^2 = R^2 - (R - halfTolerance)^2 =
439 G4double EPS = A * A * fR * kCarTolerance; /
440 if (D <= EPS) return kInfinity; // no inters
441
442 // Find intersection with Z planes
443 //
444 G4double invz = (vz == 0) ? DBL_MAX : -1./v
445 G4double dz = std::copysign(dzcut, invz);
446 G4double tzmin = (pzcut - dz) * invz;
447 G4double tzmax = (pzcut + dz) * invz;
448
449 // Find intersection with lateral surface
450 //
451 G4double tmp = -B - std::copysign(std::sqrt(
452 G4double t1 = tmp / A;
453 G4double t2 = C / tmp;
454 G4double trmin = std::min(t1, t2);
455 G4double trmax = std::max(t1, t2);
456
457 // Return distance
458 //
459 G4double tmin = std::max(tzmin, trmin);
460 G4double tmax = std::min(tzmax, trmax);
461
462 if (tmax - tmin <= halfTolerance) return kIn
463 return (tmin < halfTolerance) ? offset : tmi
464 }
465
466 //////////////////////////////////////////////
467 //
468 // Estimate distance to surface from outside
469
470 G4double G4Ellipsoid::DistanceToIn(const G4Thr
471 {
472 G4double px = p.x();
473 G4double py = p.y();
474 G4double pz = p.z();
475
476 // Safety distance to bounding box
477 G4double distX = std::abs(px) - fXmax;
478 G4double distY = std::abs(py) - fYmax;
479 G4double distZ = std::max(pz - fZTopCut, fZB
480 G4double distB = std::max(std::max(distX, di
481
482 // Safety distance to lateral surface
483 G4double x = px * fSx;
484 G4double y = py * fSy;
485 G4double z = pz * fSz;
486 G4double distR = std::sqrt(x*x + y*y + z*z)
487
488 // Return safety to in
489 G4double dist = std::max(distB, distR);
490 return (dist < 0.) ? 0. : dist;
491 }
492
493 //////////////////////////////////////////////
494 //
495 // Calculate distance to surface from inside a
496
497 G4double G4Ellipsoid::DistanceToOut(const G4Th
498 const G4Th
499 const G4bo
500 G4bo
501 G4Th
502 {
503 // Check if point flying away relative to Z
504 //
505 G4double pz = p.z() * fSz;
506 G4double vz = v.z() * fSz;
507 G4double dzcut = fZDimCut;
508 G4double pzcut = pz - fZMidCut;
509 G4double distZ = std::abs(pzcut) - dzcut;
510 if (distZ >= -halfTolerance && pzcut * vz >
511 {
512 if (calcNorm)
513 {
514 *validNorm = true;
515 n->set(0., 0., std::copysign(1., pzcut))
516 }
517 return 0.;
518 }
519
520 // Check if point is flying away relative to
521 //
522 G4double px = p.x() * fSx;
523 G4double py = p.y() * fSy;
524 G4double vx = v.x() * fSx;
525 G4double vy = v.y() * fSy;
526 G4double rr = px * px + py * py + pz * pz;
527 G4double pv = px * vx + py * vy + pz * vz;
528 G4double distR = fQ1 * rr - fQ2;
529 if (distR >= -halfTolerance && pv > 0.)
530 {
531 if (calcNorm)
532 {
533 *validNorm = true;
534 *n = G4ThreeVector(px*fSx, py*fSy, pz*fS
535 }
536 return 0.;
537 }
538
539 // Just in case check if point is outside (n
540 //
541 if (std::max(distZ, distR) > halfTolerance)
542 {
543 #ifdef G4SPECSDEBUG
544 std::ostringstream message;
545 G4long oldprc = message.precision(16);
546 message << "Point p is outside (!?) of sol
547 << GetName() << G4endl;
548 message << "Position: " << p << G4endl;;
549 message << "Direction: " << v;
550 G4cout.precision(oldprc);
551 G4Exception("G4Ellipsoid::DistanceToOut(p,
552 JustWarning, message );
553 DumpInfo();
554 #endif
555 if (calcNorm)
556 {
557 *validNorm = true;
558 *n = ApproxSurfaceNormal(p);
559 }
560 return 0.;
561 }
562
563 // Set coefficients of quadratic equation: A
564 //
565 G4double A = vx * vx + vy * vy + vz * vz;
566 G4double B = pv;
567 G4double C = rr - fR * fR;
568 G4double D = B * B - A * C;
569 // It is expected that the point is located
570 // max term in the expression for discrimina
571 // max calculation error can be derived as f
572 // A * (1 + 2e) * R^2 * (1 + 2e) = A * R^2 +
573 G4double EPS = 4. * A * fR * fR * DBL_EPSILO
574
575 if (D <= EPS) // no intersection
576 {
577 if (calcNorm)
578 {
579 *validNorm = true;
580 *n = G4ThreeVector(px*fSx, py*fSy, pz*fS
581 }
582 return 0.;
583 }
584
585 // Find intersection with Z cuts
586 //
587 G4double tzmax = (vz == 0.) ? DBL_MAX : (std
588
589 // Find intersection with lateral surface
590 //
591 G4double tmp = -B - std::copysign(std::sqrt(
592 G4double trmax = (tmp < 0.) ? C/tmp : tmp/A;
593
594 // Find distance and set normal, if required
595 //
596 G4double tmax = std::min(tzmax, trmax);
597 //if (tmax < halfTolerance) tmax = 0.;
598
599 if (calcNorm)
600 {
601 *validNorm = true;
602 if (tmax == tzmax)
603 {
604 G4double pznew = pz + tmax * vz;
605 n->set(0., 0., (pznew > fZMidCut) ? 1. :
606 }
607 else
608 {
609 G4double nx = (px + tmax * vx) * fSx;
610 G4double ny = (py + tmax * vy) * fSy;
611 G4double nz = (pz + tmax * vz) * fSz;
612 *n = G4ThreeVector(nx, ny, nz).unit();
613 }
614 }
615 return tmax;
616 }
617
618 //////////////////////////////////////////////
619 //
620 // Estimate distance to surface from inside
621
622 G4double G4Ellipsoid::DistanceToOut(const G4Th
623 {
624 // Safety distance in z direction
625 G4double distZ = std::min(fZTopCut - p.z(),
626
627 // Safety distance to lateral surface
628 G4double x = p.x() * fSx;
629 G4double y = p.y() * fSy;
630 G4double z = p.z() * fSz;
631 G4double distR = fR - std::sqrt(x*x + y*y +
632
633 // Return safety to out
634 G4double dist = std::min(distZ, distR);
635 return (dist < 0.) ? 0. : dist;
636 }
637
638 //////////////////////////////////////////////
639 //
640 // Return entity type
641
642 G4GeometryType G4Ellipsoid::GetEntityType() co
643 {
644 return {"G4Ellipsoid"};
645 }
646
647 //////////////////////////////////////////////
648 //
649 // Make a clone of the object
650
651 G4VSolid* G4Ellipsoid::Clone() const
652 {
653 return new G4Ellipsoid(*this);
654 }
655
656 //////////////////////////////////////////////
657 //
658 // Stream object contents to output stream
659
660 std::ostream& G4Ellipsoid::StreamInfo( std::os
661 {
662 G4long oldprc = os.precision(16);
663 os << "-------------------------------------
664 << " *** Dump for solid - " << GetName
665 << " =================================
666 << " Solid type: " << GetEntityType() <<
667 << " Parameters: \n"
668 << " semi-axis x: " << GetDx()/mm << "
669 << " semi-axis y: " << GetDy()/mm << "
670 << " semi-axis z: " << GetDz()/mm << "
671 << " lower cut in z: " << GetZBottomCu
672 << " upper cut in z: " << GetZTopCut()
673 << "-------------------------------------
674 os.precision(oldprc);
675 return os;
676 }
677
678 //////////////////////////////////////////////
679 //
680 // Return volume
681
682 G4double G4Ellipsoid::GetCubicVolume()
683 {
684 if (fCubicVolume == 0.)
685 {
686 G4double piAB_3 = CLHEP::pi * fDx * fDy /
687 fCubicVolume = 4. * piAB_3 * fDz;
688 if (fZBottomCut > -fDz)
689 {
690 G4double hbot = 1. + fZBottomCut / fDz;
691 fCubicVolume -= piAB_3 * hbot * hbot * (
692 }
693 if (fZTopCut < fDz)
694 {
695 G4double htop = 1. - fZTopCut / fDz;
696 fCubicVolume -= piAB_3 * htop * htop * (
697 }
698 }
699 return fCubicVolume;
700 }
701
702 //////////////////////////////////////////////
703 //
704 // Calculate area of lateral surface
705
706 G4double G4Ellipsoid::LateralSurfaceArea() con
707 {
708 constexpr G4int NPHI = 1000.;
709 constexpr G4double dPhi = CLHEP::halfpi/NPHI
710 constexpr G4double eps = 4.*DBL_EPSILON;
711
712 G4double aa = fDx*fDx;
713 G4double bb = fDy*fDy;
714 G4double cc = fDz*fDz;
715 G4double ab = fDx*fDy;
716 G4double cc_aa = cc/aa;
717 G4double cc_bb = cc/bb;
718 G4double zmax = std::min(fZTopCut, fDz);
719 G4double zmin = std::max(fZBottomCut,-fDz);
720 G4double zmax_c = zmax/fDz;
721 G4double zmin_c = zmin/fDz;
722 G4double area = 0.;
723
724 if (aa == bb) // spheroid, use analytical ex
725 {
726 G4double k = fDz/fDx;
727 G4double kk = k*k;
728 if (kk < 1. - eps)
729 {
730 G4double invk = fDx/fDz;
731 G4double root = std::sqrt(1. - kk);
732 G4double tmax = zmax_c*root;
733 G4double tmin = zmin_c*root;
734 area = CLHEP::pi*ab*
735 ((zmax_c*std::sqrt(kk + tmax*tmax) - z
736 (std::asinh(tmax*invk) - std::asinh(t
737 }
738 else if (kk > 1. + eps)
739 {
740 G4double invk = fDx/fDz;
741 G4double root = std::sqrt(kk - 1.);
742 G4double tmax = zmax_c*root;
743 G4double tmin = zmin_c*root;
744 area = CLHEP::pi*ab*
745 ((zmax_c*std::sqrt(kk - tmax*tmax) - z
746 (std::asin(tmax*invk) - std::asin(tmi
747 }
748 else
749 {
750 area = CLHEP::twopi*fDx*(zmax - zmin);
751 }
752 return area;
753 }
754
755 // ellipsoid, integration along phi
756 for (G4int i = 0; i < NPHI; ++i)
757 {
758 G4double sinPhi = std::sin(dPhi*(i + 0.5))
759 G4double kk = cc_aa + (cc_bb - cc_aa)*sinP
760 if (kk < 1. - eps)
761 {
762 G4double root = std::sqrt(1. - kk);
763 G4double tmax = zmax_c*root;
764 G4double tmin = zmin_c*root;
765 G4double invk = 1./std::sqrt(kk);
766 area += 2.*ab*dPhi*
767 ((zmax_c*std::sqrt(kk + tmax*tmax) - z
768 (std::asinh(tmax*invk) - std::asinh(t
769 }
770 else if (kk > 1. + eps)
771 {
772 G4double root = std::sqrt(kk - 1.);
773 G4double tmax = zmax_c*root;
774 G4double tmin = zmin_c*root;
775 G4double invk = 1./std::sqrt(kk);
776 area += 2.*ab*dPhi*
777 ((zmax_c*std::sqrt(kk - tmax*tmax) - z
778 (std::asin(tmax*invk) - std::asin(tmi
779 }
780 else
781 {
782 area += 4.*ab*dPhi*(zmax_c - zmin_c);
783 }
784 }
785 return area;
786 }
787
788 //////////////////////////////////////////////
789 //
790 // Return surface area
791
792 G4double G4Ellipsoid::GetSurfaceArea()
793 {
794 if (fSurfaceArea == 0.)
795 {
796 G4double piAB = CLHEP::pi * fDx * fDy;
797 fSurfaceArea = LateralSurfaceArea();
798 if (fZBottomCut > -fDz)
799 {
800 G4double hbot = 1. + fZBottomCut / fDz;
801 fSurfaceArea += piAB * hbot * (2. - hbot
802 }
803 if (fZTopCut < fDz)
804 {
805 G4double htop = 1. - fZTopCut / fDz;
806 fSurfaceArea += piAB * htop * (2. - htop
807 }
808 }
809 return fSurfaceArea;
810 }
811
812 //////////////////////////////////////////////
813 //
814 // Return random point on surface
815
816 G4ThreeVector G4Ellipsoid::GetPointOnSurface()
817 {
818 G4double A = GetDx();
819 G4double B = GetDy();
820 G4double C = GetDz();
821 G4double Zbot = GetZBottomCut();
822 G4double Ztop = GetZTopCut();
823
824 // Calculate cut areas
825 G4double Hbot = 1. + Zbot / C;
826 G4double Htop = 1. - Ztop / C;
827 G4double piAB = CLHEP::pi * A * B;
828 G4double Sbot = piAB * Hbot * (2. - Hbot);
829 G4double Stop = piAB * Htop * (2. - Htop);
830
831 // Get area of lateral surface
832 if (fLateralArea == 0.)
833 {
834 G4AutoLock l(&lateralareaMutex);
835 fLateralArea = LateralSurfaceArea();
836 l.unlock();
837 }
838 G4double Slat = fLateralArea;
839
840 // Select surface (0 - bottom cut, 1 - later
841 G4double select = (Sbot + Slat + Stop) * G4Q
842 G4int k = 0;
843 if (select > Sbot) k = 1;
844 if (select > Sbot + Slat) k = 2;
845
846 // Pick random point on selected surface (re
847 G4ThreeVector p;
848 switch (k)
849 {
850 case 0: // bootom z-cut
851 {
852 G4double scale = std::sqrt(Hbot * (2. -
853 G4TwoVector rho = G4RandomPointInEllipse
854 p.set(rho.x(), rho.y(), Zbot);
855 break;
856 }
857 case 1: // lateral surface
858 {
859 G4double x, y, z;
860 G4double mu_max = std::max(std::max(A *
861 for (G4int i = 0; i < 1000; ++i)
862 {
863 // generate random point on unit spher
864 z = (Zbot + (Ztop - Zbot) * G4QuickRan
865 G4double rho = std::sqrt((1. + z) * (1
866 G4double phi = CLHEP::twopi * G4QuickR
867 x = rho * std::cos(phi);
868 y = rho * std::sin(phi);
869 // check acceptance
870 G4double xbc = x * B * C;
871 G4double yac = y * A * C;
872 G4double zab = z * A * B;
873 G4double mu = std::sqrt(xbc * xbc + y
874 if (mu_max * G4QuickRand() <= mu) brea
875 }
876 p.set(A * x, B * y, C * z);
877 break;
878 }
879 case 2: // top z-cut
880 {
881 G4double scale = std::sqrt(Htop * (2. -
882 G4TwoVector rho = G4RandomPointInEllipse
883 p.set(rho.x(), rho.y(), Ztop);
884 break;
885 }
886 }
887 return p;
888 }
889
890 //////////////////////////////////////////////
891 //
892 // Methods for visualisation
893
894 void G4Ellipsoid::DescribeYourselfTo (G4VGraph
895 {
896 scene.AddSolid(*this);
897 }
898
899 //////////////////////////////////////////////
900 //
901 // Return vis extent
902
903 G4VisExtent G4Ellipsoid::GetExtent() const
904 {
905 return { -fXmax, fXmax, -fYmax, fYmax, fZBot
906 }
907
908 //////////////////////////////////////////////
909 //
910 // Create polyhedron
911
912 G4Polyhedron* G4Ellipsoid::CreatePolyhedron ()
913 {
914 return new G4PolyhedronEllipsoid(fDx, fDy, f
915 }
916
917 //////////////////////////////////////////////
918 //
919 // Return pointer to polyhedron
920
921 G4Polyhedron* G4Ellipsoid::GetPolyhedron () co
922 {
923 if (fpPolyhedron == nullptr ||
924 fRebuildPolyhedron ||
925 fpPolyhedron->GetNumberOfRotationStepsAt
926 fpPolyhedron->GetNumberOfRotationSteps()
927 {
928 G4AutoLock l(&polyhedronMutex);
929 delete fpPolyhedron;
930 fpPolyhedron = CreatePolyhedron();
931 fRebuildPolyhedron = false;
932 l.unlock();
933 }
934 return fpPolyhedron;
935 }
936
937 #endif // !defined(G4GEOM_USE_UELLIPSOID) || !
938