| Geant4 Cross Reference |
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 // class G4Ellipsoid
27 //
28 // Implementation of G4Ellipsoid class
29 //
30 // 10.11.99 G.Horton-Smith: first writing, based on G4Sphere class
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) && defined(G4GEOM_USE_SYS_USOLIDS))
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_INITIALIZER;
58 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZER;
59 }
60
61 using namespace CLHEP;
62
63 //////////////////////////////////////////////////////////////////////////
64 //
65 // Constructor
66
67 G4Ellipsoid::G4Ellipsoid(const G4String& name,
68 G4double xSemiAxis,
69 G4double ySemiAxis,
70 G4double zSemiAxis,
71 G4double zBottomCut,
72 G4double zTopCut)
73 : G4VSolid(name), fDx(xSemiAxis), fDy(ySemiAxis), fDz(zSemiAxis),
74 fZBottomCut(zBottomCut), fZTopCut(zTopCut)
75 {
76 CheckParameters();
77 }
78
79 //////////////////////////////////////////////////////////////////////////
80 //
81 // Fake default constructor - sets only member data and allocates memory
82 // for usage restricted to object persistency
83
84 G4Ellipsoid::G4Ellipsoid( __void__& a )
85 : G4VSolid(a), fDx(0.), fDy(0.), fDz(0.), fZBottomCut(0.), fZTopCut(0.)
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& rhs)
103 : G4VSolid(rhs),
104 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
105 fZBottomCut(rhs.fZBottomCut), fZTopCut(rhs.fZTopCut),
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.fZDimCut),
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 G4Ellipsoid& rhs)
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 tolerance
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) dimensions for Solid: "
178 << GetName() << "\n"
179 << " semi-axis x: " << fDx << "\n"
180 << " semi-axis y: " << fDy << "\n"
181 << " semi-axis z: " << fDz;
182 G4Exception("G4Ellipsoid::CheckParameters()", "GeomSolids0002",
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 || fZBottomCut >= fZTopCut)
197 {
198 std::ostringstream message;
199 message << "Invalid Z cuts for Solid: "
200 << GetName() << "\n"
201 << " bottom cut: " << fZBottomCut << "\n"
202 << " top cut: " << fZTopCut;
203 G4Exception("G4Ellipsoid::CheckParameters()", "GeomSolids0002",
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) * (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) * (1 + ratio));
224 fXmax *= scale;
225 fYmax *= scale;
226 }
227
228 // Set scale factors
229 fRsph = std::max(std::max(A, B), C); // bounding sphere
230 fR = std::min(std::min(A, B), C); // radius of sphere after scaling
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) * fSz; // middle position
237 fZDimCut = 0.5 * (fZTopCut - fZBottomCut) * fSz; // half distance
238
239 // Coefficients for approximation of distance: Q1 * (x^2 + y^2 + z^2) - Q2
240 fQ1 = 0.5 / fR;
241 fQ2 = 0.5 * fR + halfTolerance * halfTolerance * fQ1;
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 replication mechanism dimension
251 // computation & modification
252
253 void G4Ellipsoid::ComputeDimensions(G4VPVParameterisation* p,
254 const G4int n,
255 const G4VPhysicalVolume* pRep)
256 {
257 p->ComputeDimensions(*this,n,pRep);
258 }
259
260 //////////////////////////////////////////////////////////////////////////
261 //
262 // Get bounding box
263
264 void G4Ellipsoid::BoundingLimits(G4ThreeVector& pMin,
265 G4ThreeVector& pMax) const
266 {
267 pMin.set(-fXmax,-fYmax, fZBottomCut);
268 pMax.set( fXmax, fYmax, fZTopCut);
269 }
270
271 //////////////////////////////////////////////////////////////////////////
272 //
273 // Calculate extent under transform and specified limits
274
275 G4bool
276 G4Ellipsoid::CalculateExtent(const EAxis pAxis,
277 const G4VoxelLimits& pVoxelLimit,
278 const G4AffineTransform& pTransform,
279 G4double& pMin, G4double& pMax) const
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,pVoxelLimit,pTransform,pMin,pMax);
289 }
290
291 //////////////////////////////////////////////////////////////////////////
292 //
293 // Return position of point: inside/outside/on surface
294
295 EInside G4Ellipsoid::Inside(const G4ThreeVector& p) const
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) - fZDimCut;
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 : kInside;
307 }
308
309 //////////////////////////////////////////////////////////////////////////
310 //
311 // Return unit normal to surface at p
312
313 G4ThreeVector G4Ellipsoid::SurfaceNormal( const G4ThreeVector& p) const
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) - fZDimCut;
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).unit();
335 ++nsurf;
336 }
337
338 // Return normal
339 if (nsurf == 1) return norm;
340 else if (nsurf > 1) return norm.unit(); // edge
341 else
342 {
343 #ifdef G4SPECSDEBUG
344 std::ostringstream message;
345 G4long oldprc = message.precision(16);
346 message << "Point p is not on surface (!?) of solid: "
347 << GetName() << "\n";
348 message << "Position:\n";
349 message << " p.x() = " << p.x()/mm << " mm\n";
350 message << " p.y() = " << p.y()/mm << " mm\n";
351 message << " p.z() = " << p.z()/mm << " mm";
352 G4cout.precision(oldprc);
353 G4Exception("G4Ellipsoid::SurfaceNormal(p)", "GeomSolids1002",
354 JustWarning, message );
355 DumpInfo();
356 #endif
357 return ApproxSurfaceNormal(p);
358 }
359 }
360
361 //////////////////////////////////////////////////////////////////////////
362 //
363 // Find surface nearest to point and return corresponding normal.
364 // This method normally should not be called.
365
366 G4ThreeVector G4Ellipsoid::ApproxSurfaceNormal(const G4ThreeVector& p) const
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) - fZDimCut;
373 G4double distR = std::sqrt(rr) - fR;
374 if (distR > distZ && rr > 0.) // distR > distZ is correct!
375 return G4ThreeVector(x*fSx, y*fSy, z*fSz).unit();
376 else
377 return { 0., 0., std::copysign(1., z - fZMidCut) };
378 }
379
380 //////////////////////////////////////////////////////////////////////////
381 //
382 // Calculate distance to shape from outside along normalised vector
383
384 G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p,
385 const G4ThreeVector& v) const
386 {
387 G4double offset = 0.;
388 G4ThreeVector pcur = p;
389
390 // Check if point is flying away, relative to bounding box
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() >= 0.) return kInfinity;
398 if (safey >= -halfTolerance && p.y() * v.y() >= 0.) return kInfinity;
399 if (safet >= -halfTolerance && v.z() >= 0.) return kInfinity;
400 if (safeb >= -halfTolerance && v.z() <= 0.) return kInfinity;
401
402 // Relocate point, if required
403 //
404 G4double safe = std::max(std::max(std::max(safex, safey), safet), safeb);
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 : dist + offset;
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 >= 0.) return kInfinity;
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.) return kInfinity;
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 = 2 * R * halfTolerance
439 G4double EPS = A * A * fR * kCarTolerance; // discriminant at scratching
440 if (D <= EPS) return kInfinity; // no intersection or scratching
441
442 // Find intersection with Z planes
443 //
444 G4double invz = (vz == 0) ? DBL_MAX : -1./vz;
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(D), B);
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 kInfinity; // touch or no hit
463 return (tmin < halfTolerance) ? offset : tmin + offset;
464 }
465
466 //////////////////////////////////////////////////////////////////////////
467 //
468 // Estimate distance to surface from outside
469
470 G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p) const
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, fZBottomCut - pz);
480 G4double distB = std::max(std::max(distX, distY), distZ);
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) - fR;
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 along normalised vector
496
497 G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p,
498 const G4ThreeVector& v,
499 const G4bool calcNorm,
500 G4bool* validNorm,
501 G4ThreeVector* n ) const
502 {
503 // Check if point flying away relative to Z planes
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 > 0.)
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 lateral surface
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*fSz).unit();
535 }
536 return 0.;
537 }
538
539 // Just in case check if point is outside (normally it should never be)
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 solid: "
547 << GetName() << G4endl;
548 message << "Position: " << p << G4endl;;
549 message << "Direction: " << v;
550 G4cout.precision(oldprc);
551 G4Exception("G4Ellipsoid::DistanceToOut(p,v)", "GeomSolids1002",
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 t^2 + 2B t + C = 0
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 inside the sphere, so
570 // max term in the expression for discriminant is A * R^2 and
571 // max calculation error can be derived as follows:
572 // A * (1 + 2e) * R^2 * (1 + 2e) = A * R^2 + (4 * A * R^2 * e)
573 G4double EPS = 4. * A * fR * fR * DBL_EPSILON; // calculation error
574
575 if (D <= EPS) // no intersection
576 {
577 if (calcNorm)
578 {
579 *validNorm = true;
580 *n = G4ThreeVector(px*fSx, py*fSy, pz*fSz).unit();
581 }
582 return 0.;
583 }
584
585 // Find intersection with Z cuts
586 //
587 G4double tzmax = (vz == 0.) ? DBL_MAX : (std::copysign(dzcut, vz) - pzcut) / vz;
588
589 // Find intersection with lateral surface
590 //
591 G4double tmp = -B - std::copysign(std::sqrt(D), B);
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. : -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 G4ThreeVector& p) const
623 {
624 // Safety distance in z direction
625 G4double distZ = std::min(fZTopCut - p.z(), p.z() - fZBottomCut);
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 + z*z);
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() const
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::ostream& os ) const
661 {
662 G4long oldprc = os.precision(16);
663 os << "-----------------------------------------------------------\n"
664 << " *** Dump for solid - " << GetName() << " ***\n"
665 << " ===================================================\n"
666 << " Solid type: " << GetEntityType() << "\n"
667 << " Parameters: \n"
668 << " semi-axis x: " << GetDx()/mm << " mm \n"
669 << " semi-axis y: " << GetDy()/mm << " mm \n"
670 << " semi-axis z: " << GetDz()/mm << " mm \n"
671 << " lower cut in z: " << GetZBottomCut()/mm << " mm \n"
672 << " upper cut in z: " << GetZTopCut()/mm << " mm \n"
673 << "-----------------------------------------------------------\n";
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 / 3.;
687 fCubicVolume = 4. * piAB_3 * fDz;
688 if (fZBottomCut > -fDz)
689 {
690 G4double hbot = 1. + fZBottomCut / fDz;
691 fCubicVolume -= piAB_3 * hbot * hbot * (2. * fDz - fZBottomCut);
692 }
693 if (fZTopCut < fDz)
694 {
695 G4double htop = 1. - fZTopCut / fDz;
696 fCubicVolume -= piAB_3 * htop * htop * (2. * fDz + fZTopCut);
697 }
698 }
699 return fCubicVolume;
700 }
701
702 //////////////////////////////////////////////////////////////////////////
703 //
704 // Calculate area of lateral surface
705
706 G4double G4Ellipsoid::LateralSurfaceArea() const
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 expression
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) - zmin_c*std::sqrt(kk + tmin*tmin)) +
736 (std::asinh(tmax*invk) - std::asinh(tmin*invk))*kk/root);
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) - zmin_c*std::sqrt(kk - tmin*tmin)) +
746 (std::asin(tmax*invk) - std::asin(tmin*invk))*kk/root);
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)*sinPhi*sinPhi;
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) - zmin_c*std::sqrt(kk + tmin*tmin)) +
768 (std::asinh(tmax*invk) - std::asinh(tmin*invk))*kk/root);
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) - zmin_c*std::sqrt(kk - tmin*tmin)) +
778 (std::asin(tmax*invk) - std::asin(tmin*invk))*kk/root);
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() const
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 - lateral surface, 2 - top cut)
841 G4double select = (Sbot + Slat + Stop) * G4QuickRand();
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 (rejection sampling)
847 G4ThreeVector p;
848 switch (k)
849 {
850 case 0: // bootom z-cut
851 {
852 G4double scale = std::sqrt(Hbot * (2. - Hbot));
853 G4TwoVector rho = G4RandomPointInEllipse(A * scale, B * scale);
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 * B, A * C), B * C);
861 for (G4int i = 0; i < 1000; ++i)
862 {
863 // generate random point on unit sphere
864 z = (Zbot + (Ztop - Zbot) * G4QuickRand()) / C;
865 G4double rho = std::sqrt((1. + z) * (1. - z));
866 G4double phi = CLHEP::twopi * G4QuickRand();
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 + yac * yac + zab * zab);
874 if (mu_max * G4QuickRand() <= mu) break;
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. - Htop));
882 G4TwoVector rho = G4RandomPointInEllipse(A * scale, B * scale);
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 (G4VGraphicsScene& scene) const
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, fZBottomCut, fZTopCut };
906 }
907
908 //////////////////////////////////////////////////////////////////////////
909 //
910 // Create polyhedron
911
912 G4Polyhedron* G4Ellipsoid::CreatePolyhedron () const
913 {
914 return new G4PolyhedronEllipsoid(fDx, fDy, fDz, fZBottomCut, fZTopCut);
915 }
916
917 //////////////////////////////////////////////////////////////////////////
918 //
919 // Return pointer to polyhedron
920
921 G4Polyhedron* G4Ellipsoid::GetPolyhedron () const
922 {
923 if (fpPolyhedron == nullptr ||
924 fRebuildPolyhedron ||
925 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
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) || !defined(G4GEOM_USE_SYS_USOLIDS)
938