| 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 G4GeomTools implementation
27 //
28 // 10.10.2016, E.Tcherniaev: initial version.
29 // -------------------------------------------
30
31 #include "G4GeomTools.hh"
32
33 #include "geomdefs.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4GeometryTolerance.hh"
36
37 //////////////////////////////////////////////
38 //
39 // Calculate area of a triangle in 2D
40
41 G4double G4GeomTools::TriangleArea(G4double Ax
42 G4double Bx
43 G4double Cx
44 {
45 return ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax))*0
46 }
47
48 //////////////////////////////////////////////
49 //
50 // Calculate area of a triangle in 2D
51
52 G4double G4GeomTools::TriangleArea(const G4Two
53 const G4Two
54 const G4Two
55 {
56 G4double Ax = A.x(), Ay = A.y();
57 return ((B.x()-Ax)*(C.y()-Ay) - (B.y()-Ay)*(
58 }
59
60 //////////////////////////////////////////////
61 //
62 // Calculate area of a quadrilateral in 2D
63
64 G4double G4GeomTools::QuadArea(const G4TwoVect
65 const G4TwoVect
66 const G4TwoVect
67 const G4TwoVect
68 {
69 return ((C.x()-A.x())*(D.y()-B.y()) - (C.y()
70 }
71
72 //////////////////////////////////////////////
73 //
74 // Calculate area of a polygon in 2D
75
76 G4double G4GeomTools::PolygonArea(const G4TwoV
77 {
78 auto n = (G4int)p.size();
79 if (n < 3) return 0.0; // degenerate polygon
80 G4double area = p[n-1].x()*p[0].y() - p[0].x
81 for(G4int i=1; i<n; ++i)
82 {
83 area += p[i-1].x()*p[i].y() - p[i].x()*p[i
84 }
85 return area*0.5;
86 }
87
88 //////////////////////////////////////////////
89 //
90 // Point inside 2D triangle
91
92 G4bool G4GeomTools::PointInTriangle(G4double A
93 G4double B
94 G4double C
95 G4double P
96
97 {
98 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
99 {
100 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.
101 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.
102 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.
103 }
104 else
105 {
106 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.
107 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.
108 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.
109 }
110 return true;
111 }
112
113 //////////////////////////////////////////////
114 //
115 // Point inside 2D triangle
116
117 G4bool G4GeomTools::PointInTriangle(const G4Tw
118 const G4Tw
119 const G4Tw
120 const G4Tw
121 {
122 G4double Ax = A.x(), Ay = A.y();
123 G4double Bx = B.x(), By = B.y();
124 G4double Cx = C.x(), Cy = C.y();
125 G4double Px = P.x(), Py = P.y();
126 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
127 {
128 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.
129 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.
130 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.
131 }
132 else
133 {
134 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.
135 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.
136 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.
137 }
138 return true;
139 }
140
141 //////////////////////////////////////////////
142 //
143 // Point inside 2D polygon
144
145 G4bool G4GeomTools::PointInPolygon(const G4Two
146 const G4Two
147 {
148 auto Nv = (G4int)v.size();
149 G4bool in = false;
150 for (G4int i = 0, k = Nv - 1; i < Nv; k = i+
151 {
152 if ((v[i].y() > p.y()) != (v[k].y() > p.y(
153 {
154 G4double ctg = (v[k].x()-v[i].x())/(v[k]
155 in ^= static_cast<int>(p.x() < (p.y()-v[
156 }
157 }
158 return in;
159 }
160
161 //////////////////////////////////////////////
162 //
163 // Detemine whether 2D polygon is convex or no
164
165 G4bool G4GeomTools::IsConvex(const G4TwoVector
166 {
167 static const G4double kCarTolerance =
168 G4GeometryTolerance::GetInstance()->G
169
170 G4bool gotNegative = false;
171 G4bool gotPositive = false;
172 auto n = (G4int)polygon.size();
173 if (n <= 0) return false;
174 for (G4int icur=0; icur<n; ++icur)
175 {
176 G4int iprev = (icur == 0) ? n-1 : icur-1
177 G4int inext = (icur == n-1) ? 0 : icur+1
178 G4TwoVector e1 = polygon[icur] - polygon[
179 G4TwoVector e2 = polygon[inext] - polygon[
180 G4double cross = e1.x()*e2.y() - e1.y()*e2
181 if (std::abs(cross) < kCarTolerance) retur
182 if (cross < 0) gotNegative = true;
183 if (cross > 0) gotPositive = true;
184 if (gotNegative && gotPositive) return fal
185 }
186 return true;
187 }
188
189 //////////////////////////////////////////////
190 //
191 // Triangulate simple polygon
192
193 G4bool G4GeomTools::TriangulatePolygon(const G
194 G
195 {
196 result.resize(0);
197 std::vector<G4int> triangles;
198 G4bool reply = TriangulatePolygon(polygon,tr
199
200 auto n = (G4int)triangles.size();
201 for (G4int i=0; i<n; ++i) result.push_back(p
202 return reply;
203 }
204
205 //////////////////////////////////////////////
206 //
207 // Triangulation of a simple polygon by "ear c
208
209 G4bool G4GeomTools::TriangulatePolygon(const G
210 s
211 {
212 result.resize(0);
213
214 // allocate and initialize list of Vertices
215 //
216 auto n = (G4int)polygon.size();
217 if (n < 3) return false;
218
219 // we want a counter-clockwise polygon in V
220 //
221 G4double area = G4GeomTools::PolygonArea(pol
222 auto V = new G4int[n];
223 if (area > 0.)
224 for (G4int i=0; i<n; ++i) V[i] = i;
225 else
226 for (G4int i=0; i<n; ++i) V[i] = (n-1)-i;
227
228 // Triangulation: remove nv-2 Vertices, cre
229 //
230 G4int nv = n;
231 G4int count = 2*nv; // error detection count
232 for(G4int b=nv-1; nv>2; )
233 {
234 // ERROR: if we loop, it is probably a non
235 if ((count--) <= 0)
236 {
237 delete [] V;
238 if (area < 0.) std::reverse(result.begin
239 return false;
240 }
241
242 // three consecutive vertices in current p
243 G4int a = (b < nv) ? b : 0; // previou
244 b = (a+1 < nv) ? a+1 : 0; // current
245 G4int c = (b+1 < nv) ? b+1 : 0; // next
246
247 if (CheckSnip(polygon, a,b,c, nv,V))
248 {
249 // output Triangle
250 result.push_back(V[a]);
251 result.push_back(V[b]);
252 result.push_back(V[c]);
253
254 // remove vertex b from remaining polygo
255 nv--;
256 for(G4int i=b; i<nv; ++i) V[i] = V[i+1];
257
258 count = 2*nv; // resest error detection
259 }
260 }
261 delete [] V;
262 if (area < 0.) std::reverse(result.begin(),r
263 return true;
264 }
265
266 //////////////////////////////////////////////
267 //
268 // Helper function for "ear clipping" polygon
269 // Check for a valid snip
270
271 G4bool G4GeomTools::CheckSnip(const G4TwoVecto
272 G4int a, G4int b
273 G4int n, const G
274 {
275 static const G4double kCarTolerance =
276 G4GeometryTolerance::GetInstance()->G
277
278 // check orientation of Triangle
279 G4double Ax = contour[V[a]].x(), Ay = contou
280 G4double Bx = contour[V[b]].x(), By = contou
281 G4double Cx = contour[V[c]].x(), Cy = contou
282 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCar
283
284 // check that there is no point inside Trian
285 G4double xmin = std::min(std::min(Ax,Bx),Cx)
286 G4double xmax = std::max(std::max(Ax,Bx),Cx)
287 G4double ymin = std::min(std::min(Ay,By),Cy)
288 G4double ymax = std::max(std::max(Ay,By),Cy)
289 for (G4int i=0; i<n; ++i)
290 {
291 if((i == a) || (i == b) || (i == c)) conti
292 G4double Px = contour[V[i]].x();
293 if (Px < xmin || Px > xmax) continue;
294 G4double Py = contour[V[i]].y();
295 if (Py < ymin || Py > ymax) continue;
296 if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,P
297 }
298 return true;
299 }
300
301 //////////////////////////////////////////////
302 //
303 // Remove collinear and coincident points from
304
305 void G4GeomTools::RemoveRedundantVertices(G4Tw
306 std:
307 G4do
308 {
309 iout.resize(0);
310 // set tolerance squared
311 G4double delta = sqr(tolerance);
312 // set special value to mark vertices for re
313 G4double removeIt = kInfinity;
314
315 auto nv = (G4int)polygon.size();
316
317 // Main loop: check every three consecutive
318 // are collinear then mark middle point for
319 //
320 G4int icur = 0, iprev = 0, inext = 0, nout =
321 for (G4int i=0; i<nv; ++i)
322 {
323 icur = i; // index of c
324
325 for (G4int k=1; k<nv+1; ++k) // set index
326 {
327 iprev = icur - k;
328 if (iprev < 0) iprev += nv;
329 if (polygon[iprev].x() != removeIt) brea
330 }
331
332 for (G4int k=1; k<nv+1; ++k) // set index
333 {
334 inext = icur + k;
335 if (inext >= nv) inext -= nv;
336 if (polygon[inext].x() != removeIt) brea
337 }
338
339 if (iprev == inext) break; // degenerate
340
341 // Calculate parameters of triangle (iprev
342 // if triangle is too small or too narrow
343 // point for removal
344 G4TwoVector e1 = polygon[iprev] - polygon[
345 G4TwoVector e2 = polygon[inext] - polygon[
346
347 // Check length of edges, then check heigh
348 G4double leng1 = e1.mag2();
349 G4double leng2 = e2.mag2();
350 G4double leng3 = (e2-e1).mag2();
351 if (leng1 <= delta || leng2 <= delta || le
352 {
353 polygon[icur].setX(removeIt); ++nout;
354 }
355 else
356 {
357 G4double lmax = std::max(std::max(leng1,
358 G4double area = std::abs(e1.x()*e2.y()-e
359 if (area/std::sqrt(lmax) <= std::abs(tol
360 {
361 polygon[icur].setX(removeIt); ++nout;
362 }
363 }
364 }
365
366 // Remove marked points
367 //
368 icur = 0;
369 if (nv - nout < 3) // degenerate p
370 {
371 for (G4int i=0; i<nv; ++i) iout.push_back(
372 polygon.resize(0);
373 nv = 0;
374 }
375 for (G4int i=0; i<nv; ++i) // move points, i
376 {
377 if (polygon[i].x() != removeIt)
378 polygon[icur++] = polygon[i];
379 else
380 iout.push_back(i);
381 }
382 if (icur < nv) polygon.resize(icur);
383 return;
384 }
385
386 //////////////////////////////////////////////
387 //
388 // Find bounding rectangle of a disk sector
389
390 G4bool G4GeomTools::DiskExtent(G4double rmin,
391 G4double startP
392 G4TwoVector& pm
393 {
394 static const G4double kCarTolerance =
395 G4GeometryTolerance::GetInstance()->G
396
397 // check parameters
398 //
399 pmin.set(0,0);
400 pmax.set(0,0);
401 if (rmin < 0) return f
402 if (rmax <= rmin + kCarTolerance) return f
403 if (delPhi <= 0 + kCarTolerance) return f
404
405 // calculate extent
406 //
407 pmin.set(-rmax,-rmax);
408 pmax.set( rmax, rmax);
409 if (delPhi >= CLHEP::twopi) return true;
410
411 DiskExtent(rmin,rmax,
412 std::sin(startPhi),std::cos(start
413 std::sin(startPhi+delPhi),std::co
414 pmin,pmax);
415 return true;
416 }
417
418 //////////////////////////////////////////////
419 //
420 // Find bounding rectangle of a disk sector, f
421 // No check of parameters !!!
422
423 void G4GeomTools::DiskExtent(G4double rmin, G4
424 G4double sinStart
425 G4double sinEnd,
426 G4TwoVector& pmin
427 {
428 static const G4double kCarTolerance =
429 G4GeometryTolerance::GetInstance()->G
430
431 // check if 360 degrees
432 //
433 pmin.set(-rmax,-rmax);
434 pmax.set( rmax, rmax);
435
436 if (std::abs(sinEnd-sinStart) < kCarToleranc
437 std::abs(cosEnd-cosStart) < kCarToleranc
438
439 // get start and end quadrants
440 //
441 // 1 | 0
442 // ---+---
443 // 3 | 2
444 //
445 G4int icase = (cosEnd < 0) ? 1 : 0;
446 if (sinEnd < 0) icase += 2;
447 if (cosStart < 0) icase += 4;
448 if (sinStart < 0) icase += 8;
449
450 switch (icase)
451 {
452 // start quadrant 0
453 case 0: //
454 if (sinEnd < sinStart) break;
455 pmin.set(rmin*cosEnd,rmin*sinStart);
456 pmax.set(rmax*cosStart,rmax*sinEnd );
457 break;
458 case 1: //
459 pmin.set(rmax*cosEnd,std::min(rmin*sinStar
460 pmax.set(rmax*cosStart,rmax );
461 break;
462 case 2: //
463 pmin.set(-rmax,-rmax);
464 pmax.set(std::max(rmax*cosStart,rmax*cosEn
465 break;
466 case 3: //
467 pmin.set(-rmax,rmax*sinEnd);
468 pmax.set(rmax*cosStart,rmax);
469 break;
470 // start quadrant 1
471 case 4: //
472 pmin.set(-rmax,-rmax);
473 pmax.set(rmax,std::max(rmax*sinStart,rmax*
474 break;
475 case 5: //
476 if (sinEnd > sinStart) break;
477 pmin.set(rmax*cosEnd,rmin*sinEnd );
478 pmax.set(rmin*cosStart,rmax*sinStart);
479 break;
480 case 6: //
481 pmin.set(-rmax,-rmax);
482 pmax.set(rmax*cosEnd,rmax*sinStart);
483 break;
484 case 7: //
485 pmin.set(-rmax,rmax*sinEnd);
486 pmax.set(std::max(rmin*cosStart,rmin*cosEn
487 break;
488 // start quadrant 2
489 case 8: //
490 pmin.set(std::min(rmin*cosStart,rmin*cosEn
491 pmax.set(rmax,rmax*sinEnd);
492 break;
493 case 9: //
494 pmin.set(rmax*cosEnd,rmax*sinStart);
495 pmax.set(rmax,rmax);
496 break;
497 case 10: //
498 if (sinEnd < sinStart) break;
499 pmin.set(rmin*cosStart,rmax*sinStart);
500 pmax.set(rmax*cosEnd,rmin*sinEnd );
501 break;
502 case 11: //
503 pmin.set(-rmax,std::min(rmax*sinStart,rmax
504 pmax.set(rmax,rmax);
505 break;
506 // start quadrant 3
507 case 12: //
508 pmin.set(rmax*cosStart,-rmax);
509 pmax.set(rmax,rmax*sinEnd);
510 break;
511 case 13: //
512 pmin.set(std::min(rmax*cosStart,rmax*cosEn
513 pmax.set(rmax,rmax);
514 break;
515 case 14: //
516 pmin.set(rmax*cosStart,-rmax);
517 pmax.set(rmax*cosEnd,std::max(rmin*sinStar
518 break;
519 case 15: //
520 if (sinEnd > sinStart) break;
521 pmin.set(rmax*cosStart,rmax*sinEnd);
522 pmax.set(rmin*cosEnd,rmin*sinStart);
523 break;
524 }
525 return;
526 }
527
528 //////////////////////////////////////////////
529 //
530 // Compute the circumference (perimeter) of an
531
532 G4double G4GeomTools::EllipsePerimeter(G4doubl
533 {
534 G4double x = std::abs(pA);
535 G4double y = std::abs(pB);
536 G4double a = std::max(x,y);
537 G4double b = std::min(x,y);
538 G4double e = std::sqrt((1. - b/a)*(1. + b/a)
539 return 4. * a * comp_ellint_2(e);
540 }
541
542 //////////////////////////////////////////////
543 //
544 // Compute the lateral surface area of an elli
545
546 G4double G4GeomTools::EllipticConeLateralArea(
547
548
549 {
550 G4double x = std::abs(pA);
551 G4double y = std::abs(pB);
552 G4double h = std::abs(pH);
553 G4double a = std::max(x,y);
554 G4double b = std::min(x,y);
555 G4double e = std::sqrt((1. - b/a)*(1. + b/a)
556 return 2. * a * std::hypot(b,h) * comp_ellin
557 }
558
559 //////////////////////////////////////////////
560 //
561 // Compute Elliptical Integral of the Second K
562 //
563 // The algorithm is based upon Carlson B.C., "
564 // or complex elliptic integrals", Numerical A
565 // Volume 10, Issue 1, 1995 (see equations 2.3
566 //
567 // The code was adopted from C code at:
568 // http://paulbourke.net/geometry/ellipsecirc/
569
570 G4double G4GeomTools::comp_ellint_2(G4double e
571 {
572 const G4double eps = 1. / 134217728.; // 1 /
573
574 G4double a = 1.;
575 G4double b = std::sqrt((1. - e)*(1. + e));
576 if (b == 1.) return CLHEP::halfpi;
577 if (b == 0.) return 1.;
578
579 G4double x = 1.;
580 G4double y = b;
581 G4double S = 0.;
582 G4double M = 1.;
583 while (x - y > eps*y)
584 {
585 G4double tmp = (x + y) * 0.5;
586 y = std::sqrt(x*y);
587 x = tmp;
588 M += M;
589 S += M * (x - y)*(x - y);
590 }
591 return 0.5 * CLHEP::halfpi * ((a + b)*(a + b
592 }
593
594 //////////////////////////////////////////////
595 //
596 // Calcuate area of a triangle in 3D
597
598 G4ThreeVector G4GeomTools::TriangleAreaNormal(
599
600
601 {
602 return ((B-A).cross(C-A))*0.5;
603 }
604
605 //////////////////////////////////////////////
606 //
607 // Calcuate area of a quadrilateral in 3D
608
609 G4ThreeVector G4GeomTools::QuadAreaNormal(cons
610 cons
611 cons
612 cons
613 {
614 return ((C-A).cross(D-B))*0.5;
615 }
616
617 //////////////////////////////////////////////
618 //
619 // Calculate area of a polygon in 3D
620
621 G4ThreeVector G4GeomTools::PolygonAreaNormal(c
622 {
623 auto n = (G4int)p.size();
624 if (n < 3) return {0,0,0}; // degerate polyg
625 G4ThreeVector normal = p[n-1].cross(p[0]);
626 for(G4int i=1; i<n; ++i)
627 {
628 normal += p[i-1].cross(p[i]);
629 }
630 return normal*0.5;
631 }
632
633 //////////////////////////////////////////////
634 //
635 // Calculate distance between point P and line
636
637 G4double G4GeomTools::DistancePointSegment(con
638 con
639 con
640 {
641 G4ThreeVector AP = P - A;
642 G4ThreeVector AB = B - A;
643
644 G4double u = AP.dot(AB);
645 if (u <= 0) return AP.mag(); // closes
646
647 G4double len2 = AB.mag2();
648 if (u >= len2) return (B-P).mag(); // closes
649
650 return ((u/len2)*AB - AP).mag(); // distan
651 }
652
653 //////////////////////////////////////////////
654 //
655 // Find closest point on line segment in 3D
656
657 G4ThreeVector
658 G4GeomTools::ClosestPointOnSegment(const G4Thr
659 const G4Thr
660 const G4Thr
661 {
662 G4ThreeVector AP = P - A;
663 G4ThreeVector AB = B - A;
664
665 G4double u = AP.dot(AB);
666 if (u <= 0) return A; // closest point
667
668 G4double len2 = AB.mag2();
669 if (u >= len2) return B; // closest point
670
671 G4double t = u/len2;
672 return A + t*AB; // closest point
673 }
674
675 //////////////////////////////////////////////
676 //
677 // Find closest point on triangle in 3D.
678 //
679 // The implementation is based on the algorith
680 // "Geometric Tools for Computer Graphics", Ph
681 // David H Eberly, Elsevier Science (USA), 200
682 //
683 // The algorithm is also available at:
684 // http://www.geometrictools.com/Documentation
685
686 G4ThreeVector
687 G4GeomTools::ClosestPointOnTriangle(const G4Th
688 const G4Th
689 const G4Th
690 const G4Th
691 {
692 G4ThreeVector diff = A - P;
693 G4ThreeVector edge0 = B - A;
694 G4ThreeVector edge1 = C - A;
695
696 G4double a = edge0.mag2();
697 G4double b = edge0.dot(edge1);
698 G4double c = edge1.mag2();
699 G4double d = diff.dot(edge0);
700 G4double e = diff.dot(edge1);
701
702 G4double det = a*c - b*b;
703 G4double t0 = b*e - c*d;
704 G4double t1 = b*d - a*e;
705
706 /*
707 ^ t1
708 \ 2 |
709 \ |
710 \ | regions
711 \|
712 C
713 |\
714 3 | \ 1
715 | \
716 | 0 \
717 | \
718 ---- A --- B ----> t0
719 | \
720 4 | 5 \ 6
721 | \
722 */
723
724 G4int region = -1;
725 if (t0+t1 <= det)
726 region = (t0 < 0) ? ((t1 < 0) ? 4 : 3) : (
727 else
728 region = (t0 < 0) ? 2 : ((t1 < 0) ? 6 : 1)
729
730 switch (region)
731 {
732 case 0: // interior of triangle
733 {
734 G4double invDet = 1./det;
735 return A + (t0*invDet)*edge0 + (t1*invDe
736 }
737 case 1: // edge BC
738 {
739 G4double numer = c + e - b - d;
740 if (numer <= 0) return C;
741 G4double denom = a - 2*b + c;
742 return (numer >= denom) ? B : C + (numer
743 }
744 case 2: // edge AC or BC
745 {
746 G4double tmp0 = b + d;
747 G4double tmp1 = c + e;
748 if (tmp1 > tmp0)
749 {
750 G4double numer = tmp1 - tmp0;
751 G4double denom = a - 2*b + c;
752 return (numer >= denom) ? B : C + (num
753 }
754 // same: (e >= 0) ? A : ((-e >= c) ? C
755 return (tmp1 <= 0) ? C : (( e >= 0) ? A
756 }
757 case 3: // edge AC
758 return (e >= 0) ? A : ((-e >= c) ? C : A +
759
760 case 4: // edge AB or AC
761 if (d < 0) return (-d >= a) ? B : A +
762 return (e >= 0) ? A : ((-e >= c) ? C : A +
763
764 case 5: // edge AB
765 return (d >= 0) ? A : ((-d >= a) ? B : A +
766
767 case 6: // edge AB or BC
768 {
769 G4double tmp0 = b + e;
770 G4double tmp1 = a + d;
771 if (tmp1 > tmp0)
772 {
773 G4double numer = tmp1 - tmp0;
774 G4double denom = a - 2*b + c;
775 return (numer >= denom) ? C : B + (num
776 }
777 // same: (d >= 0) ? A : ((-d >= a) ? B
778 return (tmp1 <= 0) ? B : (( d >= 0) ? A
779 }
780 default: // impossible case
781 return {kInfinity,kInfinity,kInfinity};
782 }
783 }
784
785 //////////////////////////////////////////////
786 //
787 // Calculate bounding box of a spherical secto
788
789 G4bool
790 G4GeomTools::SphereExtent(G4double rmin, G4dou
791 G4double startTheta,
792 G4double startPhi, G
793 G4ThreeVector& pmin,
794 {
795 static const G4double kCarTolerance =
796 G4GeometryTolerance::GetInstance()->G
797
798 // check parameters
799 //
800 pmin.set(0,0,0);
801 pmax.set(0,0,0);
802 if (rmin < 0) return
803 if (rmax <= rmin + kCarTolerance) return
804 if (delTheta <= 0 + kCarTolerance) return
805 if (delPhi <= 0 + kCarTolerance) return
806
807 G4double stheta = startTheta;
808 G4double dtheta = delTheta;
809 if (stheta < 0 && stheta > CLHEP::pi) return
810 if (stheta + dtheta > CLHEP::pi) dtheta
811 if (dtheta <= 0 + kCarTolerance) return
812
813 // calculate extent
814 //
815 pmin.set(-rmax,-rmax,-rmax);
816 pmax.set( rmax, rmax, rmax);
817 if (dtheta >= CLHEP::pi && delPhi >= CLHEP::
818
819 G4double etheta = stheta + dtheta;
820 G4double sinStart = std::sin(stheta);
821 G4double cosStart = std::cos(stheta);
822 G4double sinEnd = std::sin(etheta);
823 G4double cosEnd = std::cos(etheta);
824
825 G4double rhomin = rmin*std::min(sinStart,sin
826 G4double rhomax = rmax;
827 if (stheta > CLHEP::halfpi) rhomax = rmax*si
828 if (etheta < CLHEP::halfpi) rhomax = rmax*si
829
830 G4TwoVector xymin,xymax;
831 DiskExtent(rhomin,rhomax,
832 std::sin(startPhi),std::cos(start
833 std::sin(startPhi+delPhi),std::co
834 xymin,xymax);
835
836 G4double zmin = std::min(rmin*cosEnd,rmax*co
837 G4double zmax = std::max(rmin*cosStart,rmax*
838 pmin.set(xymin.x(),xymin.y(),zmin);
839 pmax.set(xymax.x(),xymax.y(),zmax);
840 return true;
841 }
842