ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolyPhiFace.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PolyPhiFace.cc
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 // Implementation of G4PolyPhiFace, the face that bounds a polycone or
27 // polyhedra at its phi opening.
28 //
29 // Author: David C. Williams (davidw@scipp.ucsc.edu)
30 // --------------------------------------------------------------------
31 
32 #include "G4PolyPhiFace.hh"
33 #include "G4ClippablePolygon.hh"
34 #include "G4ReduciblePolygon.hh"
35 #include "G4AffineTransform.hh"
36 #include "G4SolidExtentList.hh"
37 #include "G4GeometryTolerance.hh"
38 
39 #include "Randomize.hh"
40 #include "G4TwoVector.hh"
41 
42 // Constructor
43 //
44 // Points r,z should be supplied in clockwise order in r,z. For example:
45 //
46 // [1]---------[2] ^ R
47 // | | |
48 // | | +--> z
49 // [0]---------[3]
50 //
52  G4double phi,
53  G4double deltaPhi,
54  G4double phiOther )
55 {
57 
58  numEdges = rz->NumVertices();
59 
60  rMin = rz->Amin();
61  rMax = rz->Amax();
62  zMin = rz->Bmin();
63  zMax = rz->Bmax();
64 
65  //
66  // Is this the "starting" phi edge of the two?
67  //
68  G4bool start = (phiOther > phi);
69 
70  //
71  // Build radial vector
72  //
73  radial = G4ThreeVector( std::cos(phi), std::sin(phi), 0.0 );
74 
75  //
76  // Build normal
77  //
78  G4double zSign = start ? 1 : -1;
79  normal = G4ThreeVector( zSign*radial.y(), -zSign*radial.x(), 0 );
80 
81  //
82  // Is allBehind?
83  //
84  allBehind = (zSign*(std::cos(phiOther)*radial.y() - std::sin(phiOther)*radial.x()) < 0);
85 
86  //
87  // Adjacent edges
88  //
89  G4double midPhi = phi + (start ? +0.5 : -0.5)*deltaPhi;
90  G4double cosMid = std::cos(midPhi),
91  sinMid = std::sin(midPhi);
92  //
93  // Allocate corners
94  //
96  //
97  // Fill them
98  //
99  G4ReduciblePolygonIterator iterRZ(rz);
100 
102  G4PolyPhiFaceVertex* helper = corners;
103 
104  iterRZ.Begin();
105  do // Loop checking, 13.08.2015, G.Cosmo
106  {
107  corn->r = iterRZ.GetA();
108  corn->z = iterRZ.GetB();
109  corn->x = corn->r*radial.x();
110  corn->y = corn->r*radial.y();
111 
112  // Add pointer on prev corner
113  //
114  if( corn == corners )
115  { corn->prev = corners+numEdges-1;}
116  else
117  { corn->prev = helper; }
118 
119  // Add pointer on next corner
120  //
121  if( corn < corners+numEdges-1 )
122  { corn->next = corn+1;}
123  else
124  { corn->next = corners; }
125 
126  helper = corn;
127  } while( ++corn, iterRZ.Next() );
128 
129  //
130  // Allocate edges
131  //
133 
134  //
135  // Fill them
136  //
137  G4double rFact = std::cos(0.5*deltaPhi);
138  G4double rFactNormalize = 1.0/std::sqrt(1.0+rFact*rFact);
139 
141  * here = corners;
142  G4PolyPhiFaceEdge* edge = edges;
143  do // Loop checking, 13.08.2015, G.Cosmo
144  {
145  G4ThreeVector sideNorm;
146 
147  edge->v0 = prev;
148  edge->v1 = here;
149 
150  G4double dr = here->r - prev->r,
151  dz = here->z - prev->z;
152 
153  edge->length = std::sqrt( dr*dr + dz*dz );
154 
155  edge->tr = dr/edge->length;
156  edge->tz = dz/edge->length;
157 
158  if ((here->r < DBL_MIN) && (prev->r < DBL_MIN))
159  {
160  //
161  // Sigh! Always exceptions!
162  // This edge runs at r==0, so its adjoing surface is not a
163  // PolyconeSide or PolyhedraSide, but the opposite PolyPhiFace.
164  //
165  G4double zSignOther = start ? -1 : 1;
166  sideNorm = G4ThreeVector( zSignOther*std::sin(phiOther),
167  -zSignOther*std::cos(phiOther), 0 );
168  }
169  else
170  {
171  sideNorm = G4ThreeVector( edge->tz*cosMid,
172  edge->tz*sinMid,
173  -edge->tr*rFact );
174  sideNorm *= rFactNormalize;
175  }
176  sideNorm += normal;
177 
178  edge->norm3D = sideNorm.unit();
179  } while( edge++, prev=here, ++here < corners+numEdges );
180 
181  //
182  // Go back and fill in corner "normals"
183  //
184  G4PolyPhiFaceEdge* prevEdge = edges+numEdges-1;
185  edge = edges;
186  do // Loop checking, 13.08.2015, G.Cosmo
187  {
188  //
189  // Calculate vertex 2D normals (on the phi surface)
190  //
191  G4double rPart = prevEdge->tr + edge->tr;
192  G4double zPart = prevEdge->tz + edge->tz;
193  G4double norm = std::sqrt( rPart*rPart + zPart*zPart );
194  G4double rNorm = +zPart/norm;
195  G4double zNorm = -rPart/norm;
196 
197  edge->v0->rNorm = rNorm;
198  edge->v0->zNorm = zNorm;
199 
200  //
201  // Calculate the 3D normals.
202  //
203  // Find the vector perpendicular to the z axis
204  // that defines the plane that contains the vertex normal
205  //
206  G4ThreeVector xyVector;
207 
208  if (edge->v0->r < DBL_MIN)
209  {
210  //
211  // This is a vertex at r==0, which is a special
212  // case. The normal we will construct lays in the
213  // plane at the center of the phi opening.
214  //
215  // We also know that rNorm < 0
216  //
217  G4double zSignOther = start ? -1 : 1;
218  G4ThreeVector normalOther( zSignOther*std::sin(phiOther),
219  -zSignOther*std::cos(phiOther), 0 );
220 
221  xyVector = - normal - normalOther;
222  }
223  else
224  {
225  //
226  // This is a vertex at r > 0. The plane
227  // is the average of the normal and the
228  // normal of the adjacent phi face
229  //
230  xyVector = G4ThreeVector( cosMid, sinMid, 0 );
231  if (rNorm < 0)
232  xyVector -= normal;
233  else
234  xyVector += normal;
235  }
236 
237  //
238  // Combine it with the r/z direction from the face
239  //
240  edge->v0->norm3D = rNorm*xyVector.unit() + G4ThreeVector( 0, 0, zNorm );
241  } while( prevEdge=edge, ++edge < edges+numEdges );
242 
243  //
244  // Build point on surface
245  //
246  G4double rAve = 0.5*(rMax-rMin),
247  zAve = 0.5*(zMax-zMin);
248  surface = G4ThreeVector( rAve*radial.x(), rAve*radial.y(), zAve );
249 }
250 
251 // Diagnose
252 //
253 // Throw an exception if something is found inconsistent with
254 // the solid.
255 //
256 // For debugging purposes only
257 //
259 {
260  G4PolyPhiFaceVertex* corner = corners;
261  do // Loop checking, 13.08.2015, G.Cosmo
262  {
263  G4ThreeVector test(corner->x, corner->y, corner->z);
264  test -= 1E-6*corner->norm3D;
265 
266  if (owner->Inside(test) != kInside)
267  G4Exception( "G4PolyPhiFace::Diagnose()", "GeomSolids0002",
268  FatalException, "Bad vertex normal found." );
269  } while( ++corner < corners+numEdges );
270 }
271 
272 // Fake default constructor - sets only member data and allocates memory
273 // for usage restricted to object persistency.
274 //
276  : numEdges(0), rMin(0.), rMax(0.), zMin(0.), zMax(0.), kCarTolerance(0.)
277 {
278 }
279 
280 
281 //
282 // Destructor
283 //
285 {
286  delete [] edges;
287  delete [] corners;
288 }
289 
290 
291 //
292 // Copy constructor
293 //
295  : G4VCSGface()
296 {
297  CopyStuff( source );
298 }
299 
300 
301 //
302 // Assignment operator
303 //
305 {
306  if (this == &source) { return *this; }
307 
308  delete [] edges;
309  delete [] corners;
310 
311  CopyStuff( source );
312 
313  return *this;
314 }
315 
316 
317 //
318 // CopyStuff (protected)
319 //
321 {
322  //
323  // The simple stuff
324  //
325  numEdges = source.numEdges;
326  normal = source.normal;
327  radial = source.radial;
328  surface = source.surface;
329  rMin = source.rMin;
330  rMax = source.rMax;
331  zMin = source.zMin;
332  zMax = source.zMax;
333  allBehind = source.allBehind;
334  triangles = nullptr;
335 
336  kCarTolerance = source.kCarTolerance;
337  fSurfaceArea = source.fSurfaceArea;
338 
339  //
340  // Corner dynamic array
341  //
344  *sourceCorn = source.corners;
345  do // Loop checking, 13.08.2015, G.Cosmo
346  {
347  *corn = *sourceCorn;
348  } while( ++sourceCorn, ++corn < corners+numEdges );
349 
350  //
351  // Edge dynamic array
352  //
354 
356  * here = corners;
357  G4PolyPhiFaceEdge* edge = edges,
358  * sourceEdge = source.edges;
359  do // Loop checking, 13.08.2015, G.Cosmo
360  {
361  *edge = *sourceEdge;
362  edge->v0 = prev;
363  edge->v1 = here;
364  } while( ++sourceEdge, ++edge, prev=here, ++here < corners+numEdges );
365 }
366 
367 // Intersect
368 //
370  const G4ThreeVector& v,
371  G4bool outgoing,
372  G4double surfTolerance,
373  G4double& distance,
374  G4double& distFromSurface,
375  G4ThreeVector& aNormal,
376  G4bool& isAllBehind )
377 {
378  G4double normSign = outgoing ? +1 : -1;
379 
380  //
381  // These don't change
382  //
383  isAllBehind = allBehind;
384  aNormal = normal;
385 
386  //
387  // Correct normal? Here we have straight sides, and can safely ignore
388  // intersections where the dot product with the normal is zero.
389  //
390  G4double dotProd = normSign*normal.dot(v);
391 
392  if (dotProd <= 0) return false;
393 
394  //
395  // Calculate distance to surface. If the side is too far
396  // behind the point, we must reject it.
397  //
398  G4ThreeVector ps = p - surface;
399  distFromSurface = -normSign*ps.dot(normal);
400 
401  if (distFromSurface < -surfTolerance) return false;
402 
403  //
404  // Calculate precise distance to intersection with the side
405  // (along the trajectory, not normal to the surface)
406  //
407  distance = distFromSurface/dotProd;
408 
409  //
410  // Calculate intersection point in r,z
411  //
412  G4ThreeVector ip = p + distance*v;
413 
414  G4double r = radial.dot(ip);
415 
416  //
417  // And is it inside the r/z extent?
418  //
419  return InsideEdgesExact( r, ip.z(), normSign, p, v );
420 }
421 
422 // Distance
423 //
425 {
426  G4double normSign = outgoing ? +1 : -1;
427  //
428  // Correct normal?
429  //
430  G4ThreeVector ps = p - surface;
431  G4double distPhi = -normSign*normal.dot(ps);
432 
433  if (distPhi < -0.5*kCarTolerance)
434  return kInfinity;
435  else if (distPhi < 0)
436  distPhi = 0.0;
437 
438  //
439  // Calculate projected point in r,z
440  //
441  G4double r = radial.dot(p);
442 
443  //
444  // Are we inside the face?
445  //
446  G4double distRZ2;
447 
448  if (InsideEdges( r, p.z(), &distRZ2, 0 ))
449  {
450  //
451  // Yup, answer is just distPhi
452  //
453  return distPhi;
454  }
455  else
456  {
457  //
458  // Nope. Penalize by distance out
459  //
460  return std::sqrt( distPhi*distPhi + distRZ2 );
461  }
462 }
463 
464 // Inside
465 //
467  G4double tolerance,
468  G4double* bestDistance )
469 {
470  //
471  // Get distance along phi, which if negative means the point
472  // is nominally inside the shape.
473  //
474  G4ThreeVector ps = p - surface;
475  G4double distPhi = normal.dot(ps);
476 
477  //
478  // Calculate projected point in r,z
479  //
480  G4double r = radial.dot(p);
481 
482  //
483  // Are we inside the face?
484  //
485  G4double distRZ2;
486  G4PolyPhiFaceVertex* base3Dnorm = nullptr;
487  G4ThreeVector* head3Dnorm = nullptr;
488 
489  if (InsideEdges( r, p.z(), &distRZ2, &base3Dnorm, &head3Dnorm ))
490  {
491  //
492  // Looks like we're inside. Distance is distance in phi.
493  //
494  *bestDistance = std::fabs(distPhi);
495 
496  //
497  // Use distPhi to decide fate
498  //
499  if (distPhi < -tolerance) return kInside;
500  if (distPhi < tolerance) return kSurface;
501  return kOutside;
502  }
503  else
504  {
505  //
506  // We're outside the extent of the face,
507  // so the distance is penalized by distance from edges in RZ
508  //
509  *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
510 
511  //
512  // Use edge normal to decide fate
513  //
514  G4ThreeVector cc( base3Dnorm->r*radial.x(),
515  base3Dnorm->r*radial.y(),
516  base3Dnorm->z );
517  cc = p - cc;
518  G4double normDist = head3Dnorm->dot(cc);
519  if ( distRZ2 > tolerance*tolerance )
520  {
521  //
522  // We're far enough away that kSurface is not possible
523  //
524  return normDist < 0 ? kInside : kOutside;
525  }
526 
527  if (normDist < -tolerance) return kInside;
528  if (normDist < tolerance) return kSurface;
529  return kOutside;
530  }
531 }
532 
533 // Normal
534 //
535 // This virtual member is simple for our planer shape,
536 // which has only one normal
537 //
539  G4double* bestDistance )
540 {
541  //
542  // Get distance along phi, which if negative means the point
543  // is nominally inside the shape.
544  //
545  G4double distPhi = normal.dot(p);
546 
547  //
548  // Calculate projected point in r,z
549  //
550  G4double r = radial.dot(p);
551 
552  //
553  // Are we inside the face?
554  //
555  G4double distRZ2;
556 
557  if (InsideEdges( r, p.z(), &distRZ2, 0 ))
558  {
559  //
560  // Yup, answer is just distPhi
561  //
562  *bestDistance = std::fabs(distPhi);
563  }
564  else
565  {
566  //
567  // Nope. Penalize by distance out
568  //
569  *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
570  }
571 
572  return normal;
573 }
574 
575 // Extent
576 //
577 // This actually isn't needed by polycone or polyhedra...
578 //
580 {
582 
583  G4PolyPhiFaceVertex* corner = corners;
584  do // Loop checking, 13.08.2015, G.Cosmo
585  {
586  G4double here = axis.x()*corner->r*radial.x()
587  + axis.y()*corner->r*radial.y()
588  + axis.z()*corner->z;
589  if (here > max) max = here;
590  } while( ++corner < corners + numEdges );
591 
592  return max;
593 }
594 
595 // CalculateExtent
596 //
597 // See notes in G4VCSGface
598 //
600  const G4VoxelLimits& voxelLimit,
602  G4SolidExtentList& extentList )
603 {
604  //
605  // Construct a (sometimes big) clippable polygon,
606  //
607  // Perform the necessary transformations while doing so
608  //
609  G4ClippablePolygon polygon;
610 
611  G4PolyPhiFaceVertex* corner = corners;
612  do // Loop checking, 13.08.2015, G.Cosmo
613  {
614  G4ThreeVector point( 0, 0, corner->z );
615  point += radial*corner->r;
616 
617  polygon.AddVertexInOrder( transform.TransformPoint( point ) );
618  } while( ++corner < corners + numEdges );
619 
620  //
621  // Clip away
622  //
623  if (polygon.PartialClip( voxelLimit, axis ))
624  {
625  //
626  // Add it to the list
627  //
628  polygon.SetNormal( transform.TransformAxis(normal) );
629  extentList.AddSurface( polygon );
630  }
631 }
632 
633 // InsideEdgesExact
634 //
635 // Decide if the point in r,z is inside the edges of our face,
636 // **but** do so consistently with other faces.
637 //
638 // This routine has functionality similar to InsideEdges, but uses
639 // an algorithm to decide if a trajectory falls inside or outside the
640 // face that uses only the trajectory p,v values and the three dimensional
641 // points representing the edges of the polygon. The objective is to plug up
642 // any leaks between touching G4PolyPhiFaces (at r==0) and any other face
643 // that uses the same convention.
644 //
645 // See: "Computational Geometry in C (Second Edition)"
646 // http://cs.smith.edu/~orourke/
647 //
649  G4double normSign,
650  const G4ThreeVector& p,
651  const G4ThreeVector& v )
652 {
653  //
654  // Quick check of extent
655  //
656  if ( (r < rMin-kCarTolerance)
657  || (r > rMax+kCarTolerance) ) return false;
658 
659  if ( (z < zMin-kCarTolerance)
660  || (z > zMax+kCarTolerance) ) return false;
661 
662  //
663  // Exact check: loop over all vertices
664  //
665  G4double qx = p.x() + v.x(),
666  qy = p.y() + v.y(),
667  qz = p.z() + v.z();
668 
669  G4int answer = 0;
670  G4PolyPhiFaceVertex *corn = corners,
671  *prev = corners+numEdges-1;
672 
673  G4double cornZ, prevZ;
674 
675  prevZ = ExactZOrder( z, qx, qy, qz, v, normSign, prev );
676  do // Loop checking, 13.08.2015, G.Cosmo
677  {
678  //
679  // Get z order of this vertex, and compare to previous vertex
680  //
681  cornZ = ExactZOrder( z, qx, qy, qz, v, normSign, corn );
682 
683  if (cornZ < 0)
684  {
685  if (prevZ < 0) continue;
686  }
687  else if (cornZ > 0)
688  {
689  if (prevZ > 0) continue;
690  }
691  else
692  {
693  //
694  // By chance, we overlap exactly (within precision) with
695  // the current vertex. Continue if the same happened previously
696  // (e.g. the previous vertex had the same z value)
697  //
698  if (prevZ == 0) continue;
699 
700  //
701  // Otherwise, to decide what to do, we need to know what is
702  // coming up next. Specifically, we need to find the next vertex
703  // with a non-zero z order.
704  //
705  // One might worry about infinite loops, but the above conditional
706  // should prevent it
707  //
708  G4PolyPhiFaceVertex *next = corn;
709  G4double nextZ;
710  do // Loop checking, 13.08.2015, G.Cosmo
711  {
712  ++next;
713  if (next == corners+numEdges) next = corners;
714 
715  nextZ = ExactZOrder( z, qx, qy, qz, v, normSign, next );
716  } while( nextZ == 0 );
717 
718  //
719  // If we won't be changing direction, go to the next vertex
720  //
721  if (nextZ*prevZ < 0) continue;
722  }
723 
724 
725  //
726  // We overlap in z with the side of the face that stretches from
727  // vertex "prev" to "corn". On which side (left or right) do
728  // we lay with respect to this segment?
729  //
730  G4ThreeVector qa( qx - prev->x, qy - prev->y, qz - prev->z ),
731  qb( qx - corn->x, qy - corn->y, qz - corn->z );
732 
733  G4double aboveOrBelow = normSign*qa.cross(qb).dot(v);
734 
735  if (aboveOrBelow > 0)
736  ++answer;
737  else if (aboveOrBelow < 0)
738  --answer;
739  else
740  {
741  //
742  // A precisely zero answer here means we exactly
743  // intersect (within roundoff) the edge of the face.
744  // Return true in this case.
745  //
746  return true;
747  }
748  } while( prevZ = cornZ, prev=corn, ++corn < corners+numEdges );
749 
750  return answer!=0;
751 }
752 
753 // InsideEdges (don't care about distance)
754 //
755 // Decide if the point in r,z is inside the edges of our face
756 //
757 // This routine can be made a zillion times quicker by implementing
758 // better code, for example:
759 //
760 // int pnpoly(int npol, float *xp, float *yp, float x, float y)
761 // {
762 // int i, j, c = 0;
763 // for (i = 0, j = npol-1; i < npol; j = i++) {
764 // if ((((yp[i]<=y) && (y<yp[j])) ||
765 // ((yp[j]<=y) && (y<yp[i]))) &&
766 // (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
767 //
768 // c = !c;
769 // }
770 // return c;
771 // }
772 //
773 // See "Point in Polyon Strategies", Eric Haines [Graphic Gems IV] pp. 24-46
774 //
775 // The algorithm below is rather unique, but is based on code needed to
776 // calculate the distance to the shape. Left here for testing...
777 //
779 {
780  //
781  // Quick check of extent
782  //
783  if ( r < rMin || r > rMax ) return false;
784  if ( z < zMin || z > zMax ) return false;
785 
786  //
787  // More thorough check
788  //
789  G4double notUsed;
790 
791  return InsideEdges( r, z, &notUsed, 0 );
792 }
793 
794 // InsideEdges (care about distance)
795 //
796 // Decide if the point in r,z is inside the edges of our face
797 //
799  G4double* bestDist2,
800  G4PolyPhiFaceVertex** base3Dnorm,
801  G4ThreeVector** head3Dnorm )
802 {
803  G4double bestDistance2 = kInfinity;
804  G4bool answer = false;
805 
806  G4PolyPhiFaceEdge* edge = edges;
807  do // Loop checking, 13.08.2015, G.Cosmo
808  {
809  G4PolyPhiFaceVertex* testMe;
810  //
811  // Get distance perpendicular to the edge
812  //
813  G4double dr = (r-edge->v0->r), dz = (z-edge->v0->z);
814 
815  G4double distOut = dr*edge->tz - dz*edge->tr;
816  G4double distance2 = distOut*distOut;
817  if (distance2 > bestDistance2) continue; // No hope!
818 
819  //
820  // Check to see if normal intersects edge within the edge's boundary
821  //
822  G4double q = dr*edge->tr + dz*edge->tz;
823 
824  //
825  // If it doesn't, penalize distance2 appropriately
826  //
827  if (q < 0)
828  {
829  distance2 += q*q;
830  testMe = edge->v0;
831  }
832  else if (q > edge->length)
833  {
834  G4double s2 = q-edge->length;
835  distance2 += s2*s2;
836  testMe = edge->v1;
837  }
838  else
839  {
840  testMe = nullptr;
841  }
842 
843  //
844  // Closest edge so far?
845  //
846  if (distance2 < bestDistance2)
847  {
848  bestDistance2 = distance2;
849  if (testMe)
850  {
851  G4double distNorm = dr*testMe->rNorm + dz*testMe->zNorm;
852  answer = (distNorm <= 0);
853  if (base3Dnorm)
854  {
855  *base3Dnorm = testMe;
856  *head3Dnorm = &testMe->norm3D;
857  }
858  }
859  else
860  {
861  answer = (distOut <= 0);
862  if (base3Dnorm)
863  {
864  *base3Dnorm = edge->v0;
865  *head3Dnorm = &edge->norm3D;
866  }
867  }
868  }
869  } while( ++edge < edges + numEdges );
870 
871  *bestDist2 = bestDistance2;
872  return answer;
873 }
874 
875 // Calculation of Surface Area of a Triangle
876 // In the same time Random Point in Triangle is given
877 //
879  G4ThreeVector p2,
880  G4ThreeVector p3,
881  G4ThreeVector* p4 )
882 {
883  G4ThreeVector v, w;
884 
885  v = p3 - p1;
886  w = p1 - p2;
887  G4double lambda1 = G4UniformRand();
888  G4double lambda2 = lambda1*G4UniformRand();
889 
890  *p4=p2 + lambda1*w + lambda2*v;
891  return 0.5*(v.cross(w)).mag();
892 }
893 
894 // Compute surface area
895 //
897 {
898  if ( fSurfaceArea==0. ) { Triangulate(); }
899  return fSurfaceArea;
900 }
901 
902 // Return random point on face
903 //
905 {
906  Triangulate();
907  return surface_point;
908 }
909 
910 //
911 // Auxiliary Functions used for Finding the PointOnFace using Triangulation
912 //
913 
914 // Calculation of 2*Area of Triangle with Sign
915 //
917  G4TwoVector b,
918  G4TwoVector c )
919 {
920  return ((b.x()-a.x())*(c.y()-a.y())-
921  (c.x()-a.x())*(b.y()-a.y()));
922 }
923 
924 // Boolean function for sign of Surface
925 //
927  G4TwoVector b,
928  G4TwoVector c )
929 {
930  return Area2(a,b,c)>0;
931 }
932 
933 // Boolean function for sign of Surface
934 //
936  G4TwoVector b,
937  G4TwoVector c )
938 {
939  return Area2(a,b,c)>=0;
940 }
941 
942 // Boolean function for sign of Surface
943 //
945  G4TwoVector b,
946  G4TwoVector c )
947 {
948  return Area2(a,b,c)==0;
949 }
950 
951 // Boolean function for finding "Proper" Intersection
952 // That means Intersection of two lines segments (a,b) and (c,d)
953 //
955  G4TwoVector b,
957 {
958  if( Collinear(a,b,c) || Collinear(a,b,d)||
959  Collinear(c,d,a) || Collinear(c,d,b) ) { return false; }
960 
961  G4bool Positive;
962  Positive = !(Left(a,b,c))^!(Left(a,b,d));
963  return Positive && (!Left(c,d,a)^!Left(c,d,b));
964 }
965 
966 // Boolean function for determining if Point c is between a and b
967 // For the tree points(a,b,c) on the same line
968 //
970 {
971  if( !Collinear(a,b,c) ) { return false; }
972 
973  if(a.x()!=b.x())
974  {
975  return ((a.x()<=c.x())&&(c.x()<=b.x()))||
976  ((a.x()>=c.x())&&(c.x()>=b.x()));
977  }
978  else
979  {
980  return ((a.y()<=c.y())&&(c.y()<=b.y()))||
981  ((a.y()>=c.y())&&(c.y()>=b.y()));
982  }
983 }
984 
985 // Boolean function for finding Intersection "Proper" or not
986 // Between two line segments (a,b) and (c,d)
987 //
989  G4TwoVector b,
991 {
992  if( IntersectProp(a,b,c,d) )
993  { return true; }
994  else if( Between(a,b,c)||
995  Between(a,b,d)||
996  Between(c,d,a)||
997  Between(c,d,b) )
998  { return true; }
999  else
1000  { return false; }
1001 }
1002 
1003 // Boolean Diagonalie help to determine
1004 // if diagonal s of segment (a,b) is convex or reflex
1005 //
1008 {
1009  G4PolyPhiFaceVertex *corner = triangles;
1010  G4PolyPhiFaceVertex *corner_next=triangles;
1011 
1012  // For each Edge (corner,corner_next)
1013  do // Loop checking, 13.08.2015, G.Cosmo
1014  {
1015  corner_next=corner->next;
1016 
1017  // Skip edges incident to a of b
1018  //
1019  if( (corner!=a)&&(corner_next!=a)
1020  &&(corner!=b)&&(corner_next!=b) )
1021  {
1022  G4TwoVector rz1,rz2,rz3,rz4;
1023  rz1 = G4TwoVector(a->r,a->z);
1024  rz2 = G4TwoVector(b->r,b->z);
1025  rz3 = G4TwoVector(corner->r,corner->z);
1026  rz4 = G4TwoVector(corner_next->r,corner_next->z);
1027  if( Intersect(rz1,rz2,rz3,rz4) ) { return false; }
1028  }
1029  corner=corner->next;
1030 
1031  } while( corner != triangles );
1032 
1033  return true;
1034 }
1035 
1036 // Boolean function that determine if b is Inside Cone (a0,a,a1)
1037 // being a the center of the Cone
1038 //
1040 {
1041  // a0,a and a1 are consecutive vertices
1042  //
1043  G4PolyPhiFaceVertex *a0,*a1;
1044  a1=a->next;
1045  a0=a->prev;
1046 
1047  G4TwoVector arz,arz0,arz1,brz;
1048  arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector(a0->r,a0->z);
1049  arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVector(b->r,b->z);
1050 
1051 
1052  if(LeftOn(arz,arz1,arz0)) // If a is convex vertex
1053  {
1054  return Left(arz,brz,arz0)&&Left(brz,arz,arz1);
1055  }
1056  else // Else a is reflex
1057  {
1058  return !( LeftOn(arz,brz,arz1)&&LeftOn(brz,arz,arz0));
1059  }
1060 }
1061 
1062 // Boolean function finding if Diagonal is possible
1063 // inside Polycone or PolyHedra
1064 //
1066 {
1067  return InCone(a,b) && InCone(b,a) && Diagonalie(a,b);
1068 }
1069 
1070 // Initialisation for Triangulisation by ear tips
1071 // For details see "Computational Geometry in C" by Joseph O'Rourke
1072 //
1074 {
1075  G4PolyPhiFaceVertex* corner = triangles;
1076  G4PolyPhiFaceVertex* c_prev, *c_next;
1077 
1078  do // Loop checking, 13.08.2015, G.Cosmo
1079  {
1080  // We need to determine three consecutive vertices
1081  //
1082  c_next=corner->next;
1083  c_prev=corner->prev;
1084 
1085  // Calculation of ears
1086  //
1087  corner->ear=Diagonal(c_prev,c_next);
1088  corner=corner->next;
1089 
1090  } while( corner!=triangles );
1091 }
1092 
1093 // Triangulisation by ear tips for Polycone or Polyhedra
1094 // For details see "Computational Geometry in C" by Joseph O'Rourke
1095 //
1097 {
1098  // The copy of Polycone is made and this copy is reordered in order to
1099  // have a list of triangles. This list is used for GetPointOnFace().
1100 
1102  triangles = tri_help;
1103  G4PolyPhiFaceVertex* triang = triangles;
1104 
1105  std::vector<G4double> areas;
1106  std::vector<G4ThreeVector> points;
1107  G4double area=0.;
1108  G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
1109  v2=triangles;
1110 
1111  // Make copy for prev/next for triang=corners
1112  //
1113  G4PolyPhiFaceVertex* helper = corners;
1114  G4PolyPhiFaceVertex* helper2 = corners;
1115  do // Loop checking, 13.08.2015, G.Cosmo
1116  {
1117  triang->r = helper->r;
1118  triang->z = helper->z;
1119  triang->x = helper->x;
1120  triang->y= helper->y;
1121 
1122  // add pointer on prev corner
1123  //
1124  if( helper==corners )
1125  { triang->prev=triangles+numEdges-1; }
1126  else
1127  { triang->prev=helper2; }
1128 
1129  // add pointer on next corner
1130  //
1131  if( helper<corners+numEdges-1 )
1132  { triang->next=triang+1; }
1133  else
1134  { triang->next=triangles; }
1135  helper2=triang;
1136  helper=helper->next;
1137  triang=triang->next;
1138 
1139  } while( helper!=corners );
1140 
1141  EarInit();
1142 
1143  G4int n=numEdges;
1144  G4int i=0;
1145  G4ThreeVector p1,p2,p3,p4;
1146  const G4int max_n_loops=numEdges*10000; // protection against infinite loop
1147 
1148  // Each step of outer loop removes one ear
1149  //
1150  while(n>3) // Loop checking, 13.08.2015, G.Cosmo
1151  { // Inner loop searches for one ear
1152  v2=triangles;
1153  do // Loop checking, 13.08.2015, G.Cosmo
1154  {
1155  if(v2->ear) // Ear found. Fill variables
1156  {
1157  // (v1,v3) is diagonal
1158  //
1159  v3=v2->next; v4=v3->next;
1160  v1=v2->prev; v0=v1->prev;
1161 
1162  // Calculate areas and points
1163 
1164  p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1165  p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
1166  p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
1167 
1168  G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1169  points.push_back(p4);
1170  areas.push_back(result1);
1171  area=area+result1;
1172 
1173  // Update earity of diagonal endpoints
1174  //
1175  v1->ear=Diagonal(v0,v3);
1176  v3->ear=Diagonal(v1,v4);
1177 
1178  // Cut off the ear v2
1179  // Has to be done for a copy and not for real PolyPhiFace
1180  //
1181  v1->next=v3;
1182  v3->prev=v1;
1183  triangles=v3; // In case the head was v2
1184  --n;
1185 
1186  break; // out of inner loop
1187  } // end if ear found
1188 
1189  v2=v2->next;
1190 
1191  } while( v2!=triangles );
1192 
1193  ++i;
1194  if(i>=max_n_loops)
1195  {
1196  G4Exception( "G4PolyPhiFace::Triangulation()",
1197  "GeomSolids0003", FatalException,
1198  "Maximum number of steps is reached for triangulation!" );
1199  }
1200  } // end outer while loop
1201 
1202  if(v2->next)
1203  {
1204  // add last triangle
1205  //
1206  v2=v2->next;
1207  p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1208  p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
1209  p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
1210  G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1211  points.push_back(p4);
1212  areas.push_back(result1);
1213  area=area+result1;
1214  }
1215 
1216  // Surface Area is stored
1217  //
1218  fSurfaceArea = area;
1219 
1220  // Second Step: choose randomly one surface
1221  //
1222  G4double chose = area*G4UniformRand();
1223 
1224  // Third Step: Get a point on choosen surface
1225  //
1226  G4double Achose1, Achose2;
1227  Achose1=0; Achose2=0.;
1228  i=0;
1229  do // Loop checking, 13.08.2015, G.Cosmo
1230  {
1231  Achose2+=areas[i];
1232  if(chose>=Achose1 && chose<Achose2)
1233  {
1235  point=points[i] ;
1237  break;
1238  }
1239  i++; Achose1=Achose2;
1240  } while( i<numEdges-2 );
1241 
1242  delete [] tri_help;
1243  tri_help = nullptr;
1244 }