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