Geant4  10.01.p02
UPolyhedraSide.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UPolyhedraSide
12 //
13 // 19.09.13 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include "UUtils.hh"
18 #include <string>
19 #include <cmath>
20 #include <sstream>
21 #include "UPolyhedraSide.hh"
22 #include "UIntersectingCone.hh"
23 
24 
25 //
26 // Constructor
27 //
28 // Values for r1,z1 and r2,z2 should be specified in clockwise
29 // order in (r,z).
30 //
32  const UPolyhedraSideRZ* tail,
33  const UPolyhedraSideRZ* head,
34  const UPolyhedraSideRZ* nextRZ,
35  int theNumSide,
36  double thePhiStart,
37  double thePhiTotal,
38  bool thePhiIsOpen,
39  bool isAllBehind)
40 {
42  fSurfaceArea = 0.;
43  fPhi.first.Set(0);
44  fPhi.second = 0.0;
45 
46  //
47  // Record values
48  //
49  r[0] = tail->r;
50  z[0] = tail->z;
51  r[1] = head->r;
52  z[1] = head->z;
53 
54  double phiTotal;
55 
56  //
57  // Set phi to our convention
58  //
59  startPhi = thePhiStart;
60  while (startPhi < 0.0) startPhi += 2 * UUtils::kPi;
61 
62  phiIsOpen = thePhiIsOpen;
63  phiTotal = (phiIsOpen) ? thePhiTotal : 2 * UUtils::kPi;
64 
65  allBehind = isAllBehind;
66 
67  //
68  // Make our intersecting cone
69  //
70  cone = new UIntersectingCone(r, z);
71 
72  //
73  // Construct side plane vector Set
74  //
75  numSide = theNumSide;
76  deltaPhi = phiTotal / theNumSide;
77  endPhi = startPhi + phiTotal;
78 
80 
82 
83  //
84  // ...this is where we start
85  //
86  double phi = startPhi;
87  UVector3 a1(r[0]*std::cos(phi), r[0]*std::sin(phi), z[0]),
88  b1(r[1]*std::cos(phi), r[1]*std::sin(phi), z[1]),
89  c1(prevRZ->r * std::cos(phi), prevRZ->r * std::sin(phi), prevRZ->z),
90  d1(nextRZ->r * std::cos(phi), nextRZ->r * std::sin(phi), nextRZ->z),
91  a2, b2, c2, d2;
92  UPolyhedraSideEdge* edge = edges;
93 
94  UPolyhedraSideVec* vec = vecs;
95  do
96  {
97  //
98  // ...this is where we are going
99  //
100  phi += deltaPhi;
101  a2 = UVector3(r[0] * std::cos(phi), r[0] * std::sin(phi), z[0]);
102  b2 = UVector3(r[1] * std::cos(phi), r[1] * std::sin(phi), z[1]);
103  c2 = UVector3(prevRZ->r * std::cos(phi), prevRZ->r * std::sin(phi), prevRZ->z);
104  d2 = UVector3(nextRZ->r * std::cos(phi), nextRZ->r * std::sin(phi), nextRZ->z);
105 
106  UVector3 tt;
107 
108  //
109  // ...build some relevant vectors.
110  // the point is to sacrifice a little memory with precalcs
111  // to gain speed
112  //
113  vec->center = 0.25 * (a1 + a2 + b1 + b2);
114 
115  tt = b2 + b1 - a2 - a1;
116  vec->surfRZ = tt.Unit();
117  if (vec == vecs) lenRZ = 0.25 * tt.Mag();
118 
119  tt = b2 - b1 + a2 - a1;
120  vec->surfPhi = tt.Unit();
121  if (vec == vecs)
122  {
123  lenPhi[0] = 0.25 * tt.Mag();
124  tt = b2 - b1;
125  lenPhi[1] = (0.5 * tt.Mag() - lenPhi[0]) / lenRZ;
126  }
127 
128  tt = vec->surfPhi.Cross(vec->surfRZ);
129  vec->normal = tt.Unit();
130 
131  //
132  // ...edge normals are the average of the normals of
133  // the two faces they connect.
134  //
135  // ...edge normals are necessary if we are to accurately
136  // decide if a point is "inside" a face. For non-convex
137  // shapes, it is absolutely necessary to know information
138  // on adjacent faces to accurate determine this.
139  //
140  // ...we don't need them for the phi edges, since that
141  // information is taken care of internally. The r/z edges,
142  // however, depend on the adjacent UPolyhedraSide.
143  //
144  UVector3 a12, adj;
145 
146  a12 = a2 - a1;
147 
148  adj = 0.5 * (c1 + c2 - a1 - a2);
149  adj = adj.Cross(a12);
150  adj = adj.Unit() + vec->normal;
151  vec->edgeNorm[0] = adj.Unit();
152 
153  a12 = b1 - b2;
154  adj = 0.5 * (d1 + d2 - b1 - b2);
155  adj = adj.Cross(a12);
156  adj = adj.Unit() + vec->normal;
157  vec->edgeNorm[1] = adj.Unit();
158 
159  //
160  // ...the corners are crucial. It is important that
161  // they are calculated consistently for adjacent
162  // UPolyhedraSides, to avoid gaps caused by roundoff.
163  //
164  vec->edges[0] = edge;
165  edge->corner[0] = a1;
166  edge->corner[1] = b1;
167  edge++;
168  vec->edges[1] = edge;
169 
170  a1 = a2;
171  b1 = b2;
172  c1 = c2;
173  d1 = d2;
174  }
175  while (++vec < vecs + numSide);
176 
177  //
178  // Clean up hanging edge
179  //
180  if (phiIsOpen)
181  {
182  edge->corner[0] = a2;
183  edge->corner[1] = b2;
184  }
185  else
186  {
187  vecs[numSide - 1].edges[1] = edges;
188  }
189 
190  //
191  // Go back and fill in remaining fields in edges
192  //
193  vec = vecs;
194  UPolyhedraSideVec* prev = vecs + numSide - 1;
195  do
196  {
197  edge = vec->edges[0]; // The edge between prev and vec
198 
199  //
200  // Okay: edge normal is average of normals of adjacent faces
201  //
202  UVector3 eNorm = vec->normal + prev->normal;
203  edge->normal = eNorm.Unit();
204 
205  //
206  // Vertex normal is average of norms of adjacent surfaces (all four)
207  // However, vec->edgeNorm is Unit vector in some direction
208  // as the sum of normals of adjacent PolyhedraSide with vec.
209  // The normalization used for this vector should be the same
210  // for vec and prev.
211  //
212  eNorm = vec->edgeNorm[0] + prev->edgeNorm[0];
213  edge->cornNorm[0] = eNorm.Unit();
214 
215  eNorm = vec->edgeNorm[1] + prev->edgeNorm[1];
216  edge->cornNorm[1] = eNorm.Unit();
217  }
218  while (prev = vec, ++vec < vecs + numSide);
219 
220  if (phiIsOpen)
221  {
222  // double rFact = std::cos(0.5*deltaPhi);
223  //
224  // If phi is open, we need to patch up normals of the
225  // first and last edges and their corresponding
226  // vertices.
227  //
228  // We use vectors that are in the plane of the
229  // face. This should be safe.
230  //
231  vec = vecs;
232 
233  UVector3 normvec = vec->edges[0]->corner[0]
234  - vec->edges[0]->corner[1];
235  normvec = normvec.Cross(vec->normal);
236  if (normvec.Dot(vec->surfPhi) > 0) normvec = -normvec;
237 
238  vec->edges[0]->normal = normvec.Unit();
239 
240  vec->edges[0]->cornNorm[0] = (vec->edges[0]->corner[0]
241  - vec->center).Unit();
242  vec->edges[0]->cornNorm[1] = (vec->edges[0]->corner[1]
243  - vec->center).Unit();
244 
245  //
246  // Repeat for ending phi
247  //
248  vec = vecs + numSide - 1;
249 
250  normvec = vec->edges[1]->corner[0] - vec->edges[1]->corner[1];
251  normvec = normvec.Cross(vec->normal);
252  if (normvec.Dot(vec->surfPhi) < 0) normvec = -normvec;
253 
254  vec->edges[1]->normal = normvec.Unit();
255 
256  vec->edges[1]->cornNorm[0] = (vec->edges[1]->corner[0]
257  - vec->center).Unit();
258  vec->edges[1]->cornNorm[1] = (vec->edges[1]->corner[1]
259  - vec->center).Unit();
260  }
261 
262  //
263  // edgeNorm is the factor one multiplies the distance along vector phi
264  // on the surface of one of our sides in order to calculate the distance
265  // from the edge. (see routine DistanceAway)
266  //
267  edgeNorm = 1.0 / std::sqrt(1.0 + lenPhi[1] * lenPhi[1]);
268 }
269 
270 
271 //
272 // Fake default constructor - sets only member data and allocates memory
273 // for usage restricted to object persistency.
274 //
276  : numSide(0), startPhi(0.), deltaPhi(0.), endPhi(0.),
277  phiIsOpen(false), allBehind(false), cone(0), vecs(0), edges(0),
278  lenRZ(0.), edgeNorm(0.), kCarTolerance(0.), fSurfaceArea(0.)
279 {
280  r[0] = r[1] = 0.;
281  z[0] = z[1] = 0.;
282  lenPhi[0] = lenPhi[1] = 0.;
283 }
284 
285 
286 //
287 // Destructor
288 //
290 {
291  delete cone;
292  delete [] vecs;
293  delete [] edges;
294 }
295 
296 
297 //
298 // Copy constructor
299 //
301  : UVCSGface()
302 {
303  CopyStuff(source);
304 }
305 
306 
307 //
308 // Assignment operator
309 //
311 {
312  if (this == &source) return *this;
313 
314  delete cone;
315  delete [] vecs;
316  delete [] edges;
317 
318  CopyStuff(source);
319 
320  return *this;
321 }
322 
323 
324 //
325 // CopyStuff
326 //
328 {
329  //
330  // The simple stuff
331  //
332  numSide = source.numSide;
333  r[0] = source.r[0];
334  r[1] = source.r[1];
335  z[0] = source.z[0];
336  z[1] = source.z[1];
337  startPhi = source.startPhi;
338  deltaPhi = source.deltaPhi;
339  endPhi = source.endPhi;
340  phiIsOpen = source.phiIsOpen;
341  allBehind = source.allBehind;
342 
343  lenRZ = source.lenRZ;
344  lenPhi[0] = source.lenPhi[0];
345  lenPhi[1] = source.lenPhi[1];
346  edgeNorm = source.edgeNorm;
347 
348  kCarTolerance = source.kCarTolerance;
349  fSurfaceArea = source.fSurfaceArea;
350 
351  cone = new UIntersectingCone(*source.cone);
352 
353  //
354  // Duplicate edges
355  //
356  int numEdges = phiIsOpen ? numSide + 1 : numSide;
357  edges = new UPolyhedraSideEdge[numEdges];
358 
359  UPolyhedraSideEdge* edge = edges,
360  *sourceEdge = source.edges;
361  do
362  {
363  *edge = *sourceEdge;
364  }
365  while (++sourceEdge, ++edge < edges + numEdges);
366 
367  //
368  // Duplicate vecs
369  //
371 
372  UPolyhedraSideVec* vec = vecs,
373  *sourceVec = source.vecs;
374  do
375  {
376  *vec = *sourceVec;
377  vec->edges[0] = edges + (sourceVec->edges[0] - source.edges);
378  vec->edges[1] = edges + (sourceVec->edges[1] - source.edges);
379  }
380  while (++sourceVec, ++vec < vecs + numSide);
381 }
382 
383 
384 //
385 // Intersect
386 //
387 // Decide if a line intersects the face.
388 //
389 // Arguments:
390 // p = (in) starting point of line segment
391 // v = (in) direction of line segment (assumed a Unit vector)
392 // A, B = (in) 2d transform variables (see note top of file)
393 // normSign = (in) desired sign for Dot product with normal (see below)
394 // surfTolerance = (in) minimum distance from the surface
395 // vecs = (in) Vector Set array
396 // distance = (out) distance to surface furfilling all requirements
397 // distFromSurface = (out) distance from the surface
398 // thisNormal = (out) normal vector of the intersecting surface
399 //
400 // Return value:
401 // true if an intersection is found. Otherwise, output parameters are
402 // undefined.
403 //
404 // Notes:
405 // * normSign: if we are "inside" the shape and only want to find out how far
406 // to leave the shape, we only want to consider intersections with surfaces in
407 // which the trajectory is leaving the shape. Since the normal vectors to the
408 // surface always point outwards from the inside, this means we want the Dot
409 // product of the trajectory direction v and the normal of the side normals[i]
410 // to be positive. Thus, we should specify normSign as +1.0. Otherwise, if
411 // we are outside and want to go in, normSign should be Set to -1.0.
412 // Don't Set normSign to zero, or you will get no intersections!
413 //
414 // * surfTolerance: see notes on argument "surfTolerance" in routine
415 // "IntersectSidePlane".
416 // ----HOWEVER---- We should *not* apply this surface tolerance if the
417 // starting point is not within phi or z of the surface. Specifically,
418 // if the starting point p angle in x/y places it on a separate side from the
419 // intersection or if the starting point p is outside the z bounds of the
420 // segment, surfTolerance must be ignored or we should *always* accept the
421 // intersection!
422 // This is simply because the sides do not have infinite extent.
423 //
424 //
426  const UVector3& v,
427  bool outgoing,
428  double surfTolerance,
429  double& distance,
430  double& distFromSurface,
431  UVector3& normal,
432  bool& isAllBehind)
433 {
434  int segment = -1;
435 
436  //
437  // Testing the intersection of individual phi faces is
438  // pretty straight forward. The simple thing therefore is to
439  // form a loop and check them all in sequence.
440  //
441  // We try to quickly decide
442  // which face would be intersected. One can make a very
443  // good guess by using the intersection with a cone.
444  // However, this is only reliable in 99% of the cases.
445  //
446  // We make a decent guess as to the one or
447  // two potential faces might get intersected, and then
448  // test them. If we have the wrong face, use the test
449  // to make a better guess.
450  //
451  // Since we might have two guesses, form a queue of
452  // potential intersecting faces. Keep an array of
453  // already tested faces to avoid doing one more than
454  // once.
455  //
456  // Result: at worst, an iterative search. On average,
457  // a little more than two tests would be required.
458  //
459 
460  if (numSide > 5)
461  {
462  //todo: maybe we could even use the second solution? how much also the second would be relevant???
463  double s1, s2;
464  int solutions = cone->LineHitsCone(p, v, s1, s2);
465  if (!solutions) return false;
466  if (solutions == 2 && s2 > 0 && (s2 < s1 || s1 < 0))
467  s1 = s2;
468 
469  segment = PhiSegment(std::atan2(p.y() + s1 * v.y(), p.x() + s1 * v.x()));
470  }
471 
472  UVector3 q = p + v;
473  UVector3 ps, delta, qa, qb, qc, qd;
474 
475  int face = -1;
476  double normSign = outgoing ? 1 : -1;
477 
478  do
479  {
480  if (face == segment) continue;
481  UPolyhedraSideVec& vec = (face == -1) ? vecs[segment] : vecs[face];
482  //
483  // Correct normal?
484  //
485  double dotProd = normSign * v.Dot(vec.normal);
486  if (dotProd <= 0) continue;
487 
488  //
489  // Is this face in front of the point along the trajectory?
490  //
491  delta = p - vec.center;
492  distFromSurface = -normSign * delta.Dot(vec.normal);
493 
494  if (distFromSurface < -surfTolerance) continue;
495 
496  //
497  // phi
498  // c -------- d ^
499  // | | |
500  // a -------- b +---> r/z
501  //
502  //
503  // Do we remain on this particular segment?
504  //
505  qc = q - vec.edges[1]->corner[0];
506  qd = q - vec.edges[1]->corner[1];
507 
508  if (normSign * qc.Cross(qd).Dot(v) < 0) continue;
509 
510  qa = q - vec.edges[0]->corner[0];
511  qb = q - vec.edges[0]->corner[1];
512 
513  if (normSign * qa.Cross(qb).Dot(v) > 0) continue;
514 
515  //
516  // We found the one and only segment we might be intersecting.
517  // Do we remain within r/z bounds?
518  //
519 
520  if (r[0] > 1 / UUtils::kInfinity && normSign * qa.Cross(qc).Dot(v) < 0) return false;
521  if (r[1] > 1 / UUtils::kInfinity && normSign * qb.Cross(qd).Dot(v) > 0) return false;
522 
523  //
524  // We allow the face to be slightly behind the trajectory
525  // (surface tolerance) only if the point p is within
526  // the vicinity of the face
527  //
528  if (distFromSurface < 0)
529  {
530  ps = p - vec.center;
531 
532  double rz = ps.Dot(vec.surfRZ);
533  if (std::fabs(rz) > lenRZ + surfTolerance) return false;
534 
535  double pp = ps.Dot(vec.surfPhi);
536  if (std::fabs(pp) > lenPhi[0] + lenPhi[1]*rz + surfTolerance) return false;
537  }
538 
539  //
540  // Intersection found. Return answer.
541  //
542  distance = distFromSurface / dotProd;
543  normal = vec.normal;
544  isAllBehind = allBehind;
545  return true;
546  }
547  while (++face < numSide);
548 
549  //
550  // Oh well. Better luck next time.
551  //
552  return false;
553 }
554 
555 
556 double UPolyhedraSide::Safety(const UVector3& p, bool outgoing)
557 {
558  double normSign = outgoing ? -1 : +1;
559 
560  //
561  // Try the closest phi segment first
562  //
563  int iPhi = ClosestPhiSegment(GetPhi(p));
564 
565  UVector3 pdotc = p - vecs[iPhi].center;
566  double normDist = pdotc.Dot(vecs[iPhi].normal);
567 
568  if (normSign * normDist > -0.5 * VUSolid::Tolerance())
569  {
570  return DistanceAway(p, vecs[iPhi], &normDist);
571  }
572 
573  //
574  // Now we have an interesting problem... do we try to find the
575  // closest facing side??
576  //
577  // Considered carefully, the answer is no. We know that if we
578  // are asking for the distance out, we are supposed to be inside,
579  // and vice versa.
580  //
581 
582  return UUtils::kInfinity;
583 }
584 
585 
586 //
587 // Inside
588 //
590  double tolerance,
591  double* bestDistance)
592 {
593  //
594  // Which phi segment is closest to this point?
595  //
596  int iPhi = ClosestPhiSegment(GetPhi(p));
597 
598  double norm;
599  //
600  // Get distance to this segment
601  //
602  *bestDistance = DistanceToOneSide(p, vecs[iPhi], &norm);
603 
604  //
605  // Use distance along normal to decide return value
606  //
607  if ((std::fabs(norm) < tolerance) && (*bestDistance < 2.0 * tolerance))
608  return VUSolid::eSurface;
609 
610  if (norm < 0) return VUSolid::eInside;
611 
612  return VUSolid::eOutside;
613 }
614 
615 
616 //
617 // Normal
618 //
620  double* bestDistance)
621 {
622  //
623  // Which phi segment is closest to this point?
624  //
625  int iPhi = ClosestPhiSegment(GetPhi(p));
626 
627  //
628  // Get distance to this segment
629  //
630  double norm;
631  *bestDistance = DistanceToOneSide(p, vecs[iPhi], &norm);
632 
633  return vecs[iPhi].normal;
634 }
635 
636 
637 //
638 // Extent
639 //
641 {
642  if (axis.Perp2() < DBL_MIN)
643  {
644  //
645  // Special case
646  //
647  return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
648  }
649 
650  int iPhi, i1, i2;
651  double best;
652  UVector3* list[4];
653 
654  //
655  // Which phi segment, if any, does the axis belong to
656  //
657  iPhi = PhiSegment(GetPhi(axis));
658 
659  if (iPhi < 0)
660  {
661  //
662  // No phi segment? Check front edge of first side and
663  // last edge of second side
664  //
665  i1 = 0;
666  i2 = numSide - 1;
667  }
668  else
669  {
670  //
671  // Check all corners of matching phi side
672  //
673  i1 = iPhi;
674  i2 = iPhi;
675  }
676 
677  list[0] = vecs[i1].edges[0]->corner;
678  list[1] = vecs[i1].edges[0]->corner + 1;
679  list[2] = vecs[i2].edges[1]->corner;
680  list[3] = vecs[i2].edges[1]->corner + 1;
681 
682  //
683  // Who's biggest?
684  //
685  best = -UUtils::kInfinity;
686  UVector3** vec = list;
687  do
688  {
689  double answer = (*vec)->Dot(axis);
690  if (answer > best) best = answer;
691  }
692  while (++vec < list + 4);
693 
694  return best;
695 }
696 
697 
698 
699 //
700 // IntersectSidePlane
701 //
702 // Decide if a line correctly intersects one side plane of our segment.
703 // It is assumed that the correct side has been chosen, and thus only
704 // the z bounds (of the entire segment) are checked.
705 //
706 // normSign - To be multiplied against normal:
707 // = +1.0 normal is unchanged
708 // = -1.0 normal is reversed (now points inward)
709 //
710 // Arguments:
711 // p - (in) Point
712 // v - (in) Direction
713 // vec - (in) Description record of the side plane
714 // normSign - (in) Sign (+/- 1) to apply to normal
715 // surfTolerance - (in) Surface tolerance (generally > 0, see below)
716 // distance - (out) Distance along v to intersection
717 // distFromSurface - (out) Distance from surface normal
718 //
719 // Notes:
720 // surfTolerance - Used to decide if a point is behind the surface,
721 // a point is allow to be -surfTolerance behind the
722 // surface (as measured along the normal), but *only*
723 // if the point is within the r/z bounds + surfTolerance
724 // of the segment.
725 //
727  const UVector3& v,
728  const UPolyhedraSideVec& vec,
729  double normSign,
730  double surfTolerance,
731  double& distance,
732  double& distFromSurface)
733 {
734  //
735  // Correct normal? Here we have straight sides, and can safely ignore
736  // intersections where the Dot product with the normal is zero.
737  //
738  double dotProd = normSign * v.Dot(vec.normal);
739 
740  if (dotProd <= 0) return false;
741 
742  //
743  // Calculate distance to surface. If the side is too far
744  // behind the point, we must reject it.
745  //
746  UVector3 delta = p - vec.center;
747  distFromSurface = -normSign * delta.Dot(vec.normal);
748 
749  if (distFromSurface < -surfTolerance) return false;
750 
751  //
752  // Calculate precise distance to intersection with the side
753  // (along the trajectory, not normal to the surface)
754  //
755  distance = distFromSurface / dotProd;
756 
757  //
758  // Do we fall off the r/z extent of the segment?
759  //
760  // Calculate this very, very carefully! Why?
761  // 1. If a RZ end is at R=0, you can't miss!
762  // 2. If you just fall off in RZ, the answer must
763  // be consistent with adjacent UPolyhedraSide faces.
764  // (2) implies that only variables used by other UPolyhedraSide
765  // faces may be used, which includes only: p, v, and the edge corners.
766  // It also means that one side is a ">" or "<", which the other
767  // must be ">=" or "<=". Fortunately, this isn't a new problem.
768  // The solution below I borrowed from Joseph O'Rourke,
769  // "Computational Geometry in C (Second Edition)"
770  // See: http://cs.smith.edu/~orourke/
771  //
772  UVector3 ic = p + distance * v - vec.center;
773  double atRZ = vec.surfRZ.Dot(ic);
774 
775  if (atRZ < 0)
776  {
777  if (r[0] == 0) return true; // Can't miss!
778 
779  if (atRZ < -lenRZ * 1.2) return false; // Forget it! Missed by a mile.
780 
781  UVector3 q = p + v;
782  UVector3 qa = q - vec.edges[0]->corner[0],
783  qb = q - vec.edges[1]->corner[0];
784  UVector3 qacb = qa.Cross(qb);
785  if (normSign * qacb.Dot(v) < 0) return false;
786 
787  if (distFromSurface < 0)
788  {
789  if (atRZ < -lenRZ - surfTolerance) return false;
790  }
791  }
792  else if (atRZ > 0)
793  {
794  if (r[1] == 0) return true; // Can't miss!
795 
796  if (atRZ > lenRZ * 1.2) return false; // Missed by a mile
797 
798  UVector3 q = p + v;
799  UVector3 qa = q - vec.edges[0]->corner[1],
800  qb = q - vec.edges[1]->corner[1];
801  UVector3 qacb = qa.Cross(qb);
802  if (normSign * qacb.Dot(v) >= 0) return false;
803 
804  if (distFromSurface < 0)
805  {
806  if (atRZ > lenRZ + surfTolerance) return false;
807  }
808  }
809 
810  return true;
811 }
812 
813 
814 //
815 // LineHitsSegments
816 //
817 // Calculate which phi segments a line intersects in three dimensions.
818 // No check is made as to whether the intersections are within the z bounds of
819 // the segment.
820 //
822  const UVector3& v,
823  int* i1, int* i2)
824 {
825  double s1, s2;
826  //
827  // First, decide if and where the line intersects the cone
828  //
829  int n = cone->LineHitsCone(p, v, s1, s2);
830 
831  if (n == 0) return 0;
832 
833  //
834  // Try first intersection.
835  //
836  *i1 = PhiSegment(std::atan2(p.y() + s1 * v.y(), p.x() + s1 * v.x()));
837  if (n == 1)
838  {
839  return (*i1 < 0) ? 0 : 1;
840  }
841 
842  //
843  // Try second intersection
844  //
845  *i2 = PhiSegment(std::atan2(p.y() + s2 * v.y(), p.x() + s2 * v.x()));
846  if (*i1 == *i2) return 0;
847 
848  if (*i1 < 0)
849  {
850  if (*i2 < 0) return 0;
851  *i1 = *i2;
852  return 1;
853  }
854 
855  if (*i2 < 0) return 1;
856 
857  return 2;
858 }
859 
860 
861 //
862 // ClosestPhiSegment
863 //
864 // Decide which phi segment is closest in phi to the point.
865 // The result is the same as PhiSegment if there is no phi opening.
866 //
868 {
869  int iPhi = PhiSegment(phi0);
870  if (iPhi >= 0) return iPhi;
871 
872  //
873  // Boogers! The points falls inside the phi segment.
874  // Look for the closest point: the start, or end
875  //
876  double phi = phi0;
877 
878  while (phi < startPhi) phi += 2 * UUtils::kPi;
879  double d1 = phi - endPhi;
880 
881  while (phi > startPhi) phi -= 2 * UUtils::kPi;
882  double d2 = startPhi - phi;
883 
884  return (d2 < d1) ? 0 : numSide - 1;
885 }
886 
887 
888 //
889 // PhiSegment
890 //
891 // Decide which phi segment an angle belongs to, counting from zero.
892 // A value of -1 indicates that the phi value is outside the shape
893 // (only possible if phiTotal < 360 degrees).
894 //
896 {
897  //
898  // How far are we from phiStart? Come up with a positive answer
899  // that is less than 2*PI
900  //
901  double phi = phi0 - startPhi;
902  while (phi < 0) phi += 2 * UUtils::kPi;
903  while (phi > 2 * UUtils::kPi) phi -= 2 * UUtils::kPi;
904 
905  //
906  // Divide
907  //
908  int answer = (int)(phi / deltaPhi);
909 
910  if (answer >= numSide)
911  {
912  if (phiIsOpen)
913  return -1; // Looks like we missed
914  else
915  answer = numSide - 1; // Probably just roundoff
916  }
917 
918  return answer;
919 }
920 
921 
922 //
923 // GetPhi
924 //
925 // Calculate Phi for a given 3-vector (point), if not already cached for the
926 // same point, in the attempt to avoid consecutive computation of the same
927 // quantity
928 //
930 {
931  double val = 0.;
932 
933  if (fPhi.first != p)
934  {
935  val = p.Phi();
936  fPhi.first = p;
937  fPhi.second = val;
938  }
939  else
940  {
941  val = fPhi.second;
942  }
943  return val;
944 }
945 
946 
947 //
948 // DistanceToOneSide
949 //
950 // Arguments:
951 // p - (in) Point to check
952 // vec - (in) vector Set of this side
953 // normDist - (out) distance normal to the side or edge, as appropriate, signed
954 // Return value = total distance from the side
955 //
957  const UPolyhedraSideVec& vec,
958  double* normDist)
959 {
960  UVector3 pct = p - vec.center;
961 
962  //
963  // Get normal distance
964  //
965  *normDist = vec.normal.Dot(pct);
966 
967  //
968  // Add edge penalty
969  //
970  return DistanceAway(p, vec, normDist);
971 }
972 
973 
974 //
975 // DistanceAway
976 //
977 // Add distance from side edges, if necesssary, to total distance,
978 // and updates normDist appropriate depending on edge normals.
979 //
981  const UPolyhedraSideVec& vec,
982  double* normDist)
983 {
984  double distOut2;
985  UVector3 pct = p - vec.center;
986  double distFaceNorm = *normDist;
987  //
988  // Okay, are we inside bounds?
989  //
990  double pcDotRZ = pct.Dot(vec.surfRZ);
991  double pcDotPhi = pct.Dot(vec.surfPhi);
992 
993  //
994  // Go through all permutations.
995  // Phi
996  // | | ^
997  // B | H | E |
998  // ------[1]------------[3]----- |
999  // |XXXXXXXXXXXXXX| +----> RZ
1000  // C |XXXXXXXXXXXXXX| F
1001  // |XXXXXXXXXXXXXX|
1002  // ------[0]------------[2]----
1003  // A | G | D
1004  // | |
1005  //
1006  // It's real messy, but at least it's quick
1007  //
1008 
1009  if (pcDotRZ < -lenRZ)
1010  {
1011  double lenPhiZ = lenPhi[0] - lenRZ * lenPhi[1];
1012  double distOutZ = pcDotRZ + lenRZ;
1013  //
1014  // Below in RZ
1015  //
1016  if (pcDotPhi < -lenPhiZ)
1017  {
1018  //
1019  // ...and below in phi. Find distance to point (A)
1020  //
1021  double distOutPhi = pcDotPhi + lenPhiZ;
1022  distOut2 = distOutPhi * distOutPhi + distOutZ * distOutZ;
1023  UVector3 pa = p - vec.edges[0]->corner[0];
1024  *normDist = pa.Dot(vec.edges[0]->cornNorm[0]);
1025  }
1026  else if (pcDotPhi > lenPhiZ)
1027  {
1028  //
1029  // ...and above in phi. Find distance to point (B)
1030  //
1031  double distOutPhi = pcDotPhi - lenPhiZ;
1032  distOut2 = distOutPhi * distOutPhi + distOutZ * distOutZ;
1033  UVector3 pb = p - vec.edges[1]->corner[0];
1034  *normDist = pb.Dot(vec.edges[1]->cornNorm[0]);
1035  }
1036  else
1037  {
1038  //
1039  // ...and inside in phi. Find distance to line (C)
1040  //
1041  UVector3 pa = p - vec.edges[0]->corner[0];
1042  distOut2 = distOutZ * distOutZ;
1043  *normDist = pa.Dot(vec.edgeNorm[0]);
1044  }
1045  }
1046  else if (pcDotRZ > lenRZ)
1047  {
1048  double lenPhiZ = lenPhi[0] + lenRZ * lenPhi[1];
1049  double distOutZ = pcDotRZ - lenRZ;
1050  //
1051  // Above in RZ
1052  //
1053  if (pcDotPhi < -lenPhiZ)
1054  {
1055  //
1056  // ...and below in phi. Find distance to point (D)
1057  //
1058  double distOutPhi = pcDotPhi + lenPhiZ;
1059  distOut2 = distOutPhi * distOutPhi + distOutZ * distOutZ;
1060  UVector3 pd = p - vec.edges[0]->corner[1];
1061  *normDist = pd.Dot(vec.edges[0]->cornNorm[1]);
1062  }
1063  else if (pcDotPhi > lenPhiZ)
1064  {
1065  //
1066  // ...and above in phi. Find distance to point (E)
1067  //
1068  double distOutPhi = pcDotPhi - lenPhiZ;
1069  distOut2 = distOutPhi * distOutPhi + distOutZ * distOutZ;
1070  UVector3 pe = p - vec.edges[1]->corner[1];
1071  *normDist = pe.Dot(vec.edges[1]->cornNorm[1]);
1072  }
1073  else
1074  {
1075  //
1076  // ...and inside in phi. Find distance to line (F)
1077  //
1078  distOut2 = distOutZ * distOutZ;
1079  UVector3 pd = p - vec.edges[0]->corner[1];
1080  *normDist = pd.Dot(vec.edgeNorm[1]);
1081  }
1082  }
1083  else
1084  {
1085  double lenPhiZ = lenPhi[0] + pcDotRZ * lenPhi[1];
1086  //
1087  // We are inside RZ bounds
1088  //
1089  if (pcDotPhi < -lenPhiZ)
1090  {
1091  //
1092  // ...and below in phi. Find distance to line (G)
1093  //
1094  double distOut = edgeNorm * (pcDotPhi + lenPhiZ);
1095  distOut2 = distOut * distOut;
1096  UVector3 pd = p - vec.edges[0]->corner[1];
1097  *normDist = pd.Dot(vec.edges[0]->normal);
1098  }
1099  else if (pcDotPhi > lenPhiZ)
1100  {
1101  //
1102  // ...and above in phi. Find distance to line (H)
1103  //
1104  double distOut = edgeNorm * (pcDotPhi - lenPhiZ);
1105  distOut2 = distOut * distOut;
1106  UVector3 pe = p - vec.edges[1]->corner[1];
1107  *normDist = pe.Dot(vec.edges[1]->normal);
1108  }
1109  else
1110  {
1111  //
1112  // Inside bounds! No penalty.
1113  //
1114  return std::fabs(distFaceNorm);
1115  }
1116  }
1117  return std::sqrt(distFaceNorm * distFaceNorm + distOut2);
1118 }
1119 
1120 
1121 //
1122 // Calculation of surface area of a triangle.
1123 // At the same time a random point in the triangle is given
1124 //
1126  UVector3 p2,
1127  UVector3 p3,
1128  UVector3* p4)
1129 {
1130  UVector3 v, w;
1131 
1132  v = p3 - p1;
1133  w = p1 - p2;
1134  double lambda1 = UUtils::Random();
1135  double lambda2 = lambda1 * UUtils::Random();
1136 
1137  *p4 = p2 + lambda1 * w + lambda2 * v;
1138  return 0.5 * (v.Cross(w)).Mag();
1139 }
1140 
1141 
1142 //
1143 // GetPointOnPlane
1144 //
1145 // Auxiliary method for GetPointOnSurface()
1146 //
1147 UVector3
1149  UVector3 p2, UVector3 p3,
1150  double* Area)
1151 {
1152  double chose, aOne, aTwo;
1153  UVector3 point1, point2;
1154 
1155  aOne = SurfaceTriangle(p0, p1, p2, &point1);
1156  aTwo = SurfaceTriangle(p2, p3, p0, &point2);
1157  *Area = aOne + aTwo;
1158 
1159  chose = UUtils::Random() * (aOne + aTwo);
1160  if ((chose >= 0.) && (chose < aOne))
1161  {
1162  return (point1);
1163  }
1164  return (point2);
1165 }
1166 
1167 
1168 //
1169 // SurfaceArea()
1170 //
1172 {
1173  if (fSurfaceArea == 0.)
1174  {
1175  // Define the variables
1176  //
1177  double area, areas;
1178  UVector3 point1;
1179  UVector3 v1, v2, v3, v4;
1180  UPolyhedraSideVec* vec = vecs;
1181  areas = 0.;
1182 
1183  // Do a loop on all SideEdge
1184  //
1185  do
1186  {
1187  // Define 4points for a Plane or Triangle
1188  //
1189  v1 = vec->edges[0]->corner[0];
1190  v2 = vec->edges[0]->corner[1];
1191  v3 = vec->edges[1]->corner[1];
1192  v4 = vec->edges[1]->corner[0];
1193  point1 = GetPointOnPlane(v1, v2, v3, v4, &area);
1194  areas += area;
1195  }
1196  while (++vec < vecs + numSide);
1197 
1198  fSurfaceArea = areas;
1199  }
1200  return fSurfaceArea;
1201 }
1202 
1203 
1204 //
1205 // GetPointOnFace()
1206 //
1208 {
1209  // Define the variables
1210  //
1211  std::vector<double> areas;
1212  std::vector<UVector3> points;
1213  double area = 0;
1214  double result1;
1215  UVector3 point1;
1216  UVector3 v1, v2, v3, v4;
1217  UPolyhedraSideVec* vec = vecs;
1218 
1219  // Do a loop on all SideEdge
1220  //
1221  do
1222  {
1223  // Define 4points for a Plane or Triangle
1224  //
1225  v1 = vec->edges[0]->corner[0];
1226  v2 = vec->edges[0]->corner[1];
1227  v3 = vec->edges[1]->corner[1];
1228  v4 = vec->edges[1]->corner[0];
1229  point1 = GetPointOnPlane(v1, v2, v3, v4, &result1);
1230  points.push_back(point1);
1231  areas.push_back(result1);
1232  area += result1;
1233  }
1234  while (++vec < vecs + numSide);
1235 
1236  // Choose randomly one of the surfaces and point on it
1237  //
1238  double chose = area * UUtils::Random();
1239  double Achose1, Achose2;
1240  Achose1 = 0;
1241  Achose2 = 0.;
1242  int i = 0;
1243  do
1244  {
1245  Achose2 += areas[i];
1246  if (chose >= Achose1 && chose < Achose2)
1247  {
1248  point1 = points[i] ;
1249  break;
1250  }
1251  i++;
1252  Achose1 = Achose2;
1253  }
1254  while (i < numSide);
1255 
1256  return point1;
1257 }
bool IntersectSidePlane(const UVector3 &p, const UVector3 &v, const UPolyhedraSideVec &vec, double normSign, double surfTolerance, double &distance, double &distFromSurface)
static c2_factory< G4double > c2
double & z()
Definition: UVector3.hh:352
UIntersectingCone * cone
static const G4double d1
double & y()
Definition: UVector3.hh:348
virtual ~UPolyhedraSide()
double Phi() const
Definition: UVector3.cc:64
UVector3 Cross(const UVector3 &) const
Definition: UVector3.hh:281
static const G4double a1
UPolyhedraSideVec * vecs
double Safety(const UVector3 &p, bool outgoing)
static const G4double tolerance
double GetPhi(const UVector3 &p)
UVector3 GetPointOnPlane(UVector3 p0, UVector3 p1, UVector3 p2, UVector3 p3, double *Area)
static double Tolerance()
Definition: VUSolid.hh:127
UPolyhedraSide(const UPolyhedraSideRZ *prevRZ, const UPolyhedraSideRZ *tail, const UPolyhedraSideRZ *head, const UPolyhedraSideRZ *nextRZ, int numSide, double phiStart, double phiTotal, bool phiIsOpen, bool isAllBehind=false)
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
double SurfaceTriangle(UVector3 p1, UVector3 p2, UVector3 p3, UVector3 *p4)
UVector3 Normal(const UVector3 &p, double *bestDistance)
VUSolid::EnumInside Inside(const UVector3 &p, double tolerance, double *bestDistance)
bool Distance(const UVector3 &p, const UVector3 &v, bool outgoing, double surfTolerance, double &distance, double &distFromSurface, UVector3 &normal, bool &allBehind)
int ClosestPhiSegment(double phi)
static const double kInfinity
Definition: UUtils.hh:54
double & x()
Definition: UVector3.hh:344
const G4double p2
const G4double p1
UPolyhedraSideEdge * edges
static const G4double b2
EnumInside
Definition: VUSolid.hh:23
const G4int n
double Dot(const UVector3 &) const
Definition: UVector3.hh:276
static const G4double c1
double DistanceToOneSide(const UVector3 &p, const UPolyhedraSideVec &vec, double *normDist)
double Extent(const UVector3 axis)
double Mag() const
Definition: UVector3.cc:48
const G4double p0
double Perp2() const
Definition: UVector3.hh:292
UVector3 GetPointOnFace()
static const double kPi
Definition: UUtils.hh:49
double ZHi() const
UVector3 Unit() const
Definition: UVector3.cc:80
#define DBL_MIN
Definition: templates.hh:75
int LineHitsSegments(const UVector3 &p, const UVector3 &v, int *i1, int *i2)
static const G4double b1
static const G4double d2
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double ZLo() const
double DistanceAway(const UVector3 &p, const UPolyhedraSideVec &vec, double *normDist)
int LineHitsCone(const UVector3 &p, const UVector3 &v, double &s1, double &s2)
void CopyStuff(const UPolyhedraSide &source)
static const G4double a2
UPolyhedraSide & operator=(const UPolyhedraSide &source)
std::pair< UVector3, double > fPhi
int PhiSegment(double phi)