Geant4  10.00.p03
UPolyconeSide.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 // UPolyconeSide
12 //
13 // 19.04.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 "UPolyconeSide.hh"
22 #include "UIntersectingCone.hh"
23 #include "VUSolid.hh"
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 UPolyconeSideRZ* tail,
33  const UPolyconeSideRZ* head,
34  const UPolyconeSideRZ* nextRZ,
35  double thePhiStart,
36  double theDeltaPhi,
37  bool thePhiIsOpen,
38  bool isAllBehind)
39  : ncorners(0), corners(0)
40 {
41 
43 
44  fSurfaceArea = 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  phiIsOpen = thePhiIsOpen;
55  if (phiIsOpen)
56  {
57  deltaPhi = theDeltaPhi;
58  startPhi = thePhiStart;
59 
60  //
61  // Set phi values to our conventions
62  //
63  while (deltaPhi < 0.0) deltaPhi += 2 * UUtils::kPi;
64  while (startPhi < 0.0) startPhi += 2 * UUtils::kPi;
65 
66  //
67  // Calculate corner coordinates
68  //
69  ncorners = 4;
70  corners = new UVector3[ncorners];
71 
72  corners[0] = UVector3(tail->r * std::cos(startPhi),
73  tail->r * std::sin(startPhi), tail->z);
74  corners[1] = UVector3(head->r * std::cos(startPhi),
75  head->r * std::sin(startPhi), head->z);
76  corners[2] = UVector3(tail->r * std::cos(startPhi + deltaPhi),
77  tail->r * std::sin(startPhi + deltaPhi), tail->z);
78  corners[3] = UVector3(head->r * std::cos(startPhi + deltaPhi),
79  head->r * std::sin(startPhi + deltaPhi), head->z);
80  }
81  else
82  {
83  deltaPhi = 2 * UUtils::kPi;
84  startPhi = 0.0;
85  }
86 
87  allBehind = isAllBehind;
88 
89  //
90  // Make our intersecting cone
91  //
92  cone = new UIntersectingCone(r, z);
93 
94  //
95  // Calculate vectors in r,z space
96  //
97  rS = r[1] - r[0];
98  zS = z[1] - z[0];
99  length = std::sqrt(rS * rS + zS * zS);
100  rS /= length;
101  zS /= length;
102 
103  rNorm = +zS;
104  zNorm = -rS;
105 
106  double lAdj;
107 
108  prevRS = r[0] - prevRZ->r;
109  prevZS = z[0] - prevRZ->z;
110  lAdj = std::sqrt(prevRS * prevRS + prevZS * prevZS);
111  prevRS /= lAdj;
112  prevZS /= lAdj;
113 
114  rNormEdge[0] = rNorm + prevZS;
115  zNormEdge[0] = zNorm - prevRS;
116  lAdj = std::sqrt(rNormEdge[0] * rNormEdge[0] + zNormEdge[0] * zNormEdge[0]);
117  rNormEdge[0] /= lAdj;
118  zNormEdge[0] /= lAdj;
119 
120  nextRS = nextRZ->r - r[1];
121  nextZS = nextRZ->z - z[1];
122  lAdj = std::sqrt(nextRS * nextRS + nextZS * nextZS);
123  nextRS /= lAdj;
124  nextZS /= lAdj;
125 
126  rNormEdge[1] = rNorm + nextZS;
127  zNormEdge[1] = zNorm - nextRS;
128  lAdj = std::sqrt(rNormEdge[1] * rNormEdge[1] + zNormEdge[1] * zNormEdge[1]);
129  rNormEdge[1] /= lAdj;
130  zNormEdge[1] /= lAdj;
131 }
132 
133 //
134 // Fake default constructor - sets only member data and allocates memory
135 // for usage restricted to object persistency.
136 //
138  : startPhi(0.), deltaPhi(0.), phiIsOpen(false), allBehind(false),
139  cone(0), rNorm(0.), zNorm(0.), rS(0.), zS(0.), length(0.),
140  prevRS(0.), prevZS(0.), nextRS(0.), nextZS(0.), ncorners(0), corners(0),
141  tolerance(0.), fSurfaceArea(0.)
142 {
143  r[0] = r[1] = 0.;
144  z[0] = z[1] = 0.;
145  rNormEdge[0] = rNormEdge[1] = 0.;
146  zNormEdge[0] = zNormEdge[1] = 0.;
147 }
148 
149 //
150 // Destructor
151 //
153 {
154  delete cone;
155  if (phiIsOpen)
156  {
157  delete [] corners;
158  }
159 }
160 
161 
162 //
163 // Copy constructor
164 //
166  : UVCSGface(), ncorners(0), corners(0)
167 {
168 
169  CopyStuff(source);
170 }
171 
172 
173 //
174 // Assignment operator
175 //
177 {
178  if (this == &source)
179  {
180  return *this;
181  }
182 
183  delete cone;
184  if (phiIsOpen)
185  {
186  delete [] corners;
187  }
188 
189  CopyStuff(source);
190 
191  return *this;
192 }
193 
194 
195 //
196 // CopyStuff
197 //
199 {
200  r[0] = source.r[0];
201  r[1] = source.r[1];
202  z[0] = source.z[0];
203  z[1] = source.z[1];
204 
205  startPhi = source.startPhi;
206  deltaPhi = source.deltaPhi;
207  phiIsOpen = source.phiIsOpen;
208  allBehind = source.allBehind;
209 
210  tolerance = source.tolerance;
211  fSurfaceArea = source.fSurfaceArea;
212 
213  cone = new UIntersectingCone(*source.cone);
214 
215  rNorm = source.rNorm;
216  zNorm = source.zNorm;
217  rS = source.rS;
218  zS = source.zS;
219  length = source.length;
220  prevRS = source.prevRS;
221  prevZS = source.prevZS;
222  nextRS = source.nextRS;
223  nextZS = source.nextZS;
224 
225  rNormEdge[0] = source.rNormEdge[0];
226  rNormEdge[1] = source.rNormEdge[1];
227  zNormEdge[0] = source.zNormEdge[0];
228  zNormEdge[1] = source.zNormEdge[1];
229 
230  if (phiIsOpen)
231  {
232  ncorners = 4;
233  corners = new UVector3[ncorners];
234 
235  corners[0] = source.corners[0];
236  corners[1] = source.corners[1];
237  corners[2] = source.corners[2];
238  corners[3] = source.corners[3];
239  }
240 }
241 
242 
243 //
244 // Intersect
245 //
247  const UVector3& v,
248  bool outgoing,
249  double surfTolerance,
250  double& distance,
251  double& distFromSurface,
252  UVector3& normal,
253  bool& isAllBehind)
254 {
255  double s1, s2;
256  double normSign = outgoing ? +1 : -1;
257 
258  isAllBehind = allBehind;
259 
260  //
261  // Check for two possible intersections
262  //
263  int nside = cone->LineHitsCone(p, v, s1, s2);
264  if (nside == 0) return false;
265 
266  //
267  // Check the first side first, since it is (supposed to be) closest
268  //
269  UVector3 hit = p + s1 * v;
270 
271  if (PointOnCone(hit, normSign, p, v, normal))
272  {
273  //
274  // Good intersection! What about the normal?
275  //
276  if (normSign * v.Dot(normal) > 0)
277  {
278  //
279  // We have a valid intersection, but it could very easily
280  // be behind the point. To decide if we tolerate this,
281  // we have to see if the point p is on the surface near
282  // the intersecting point.
283  //
284  // What does it mean exactly for the point p to be "near"
285  // the intersection? It means that if we draw a line from
286  // p to the hit, the line remains entirely within the
287  // tolerance bounds of the cone. To test this, we can
288  // ask if the normal is correct near p.
289  //
290  double pr = p.Perp();
291  if (pr < DBL_MIN) pr = DBL_MIN;
292  UVector3 pNormal(rNorm * p.x / pr, rNorm * p.y / pr, zNorm);
293  if (normSign * v.Dot(pNormal) > 0)
294  {
295  //
296  // p and intersection in same hemisphere
297  //
298  double distOutside2;
299  distFromSurface = -normSign * DistanceAway(p, false, distOutside2);
300  if (distOutside2 < surfTolerance * surfTolerance)
301  {
302  if (distFromSurface > -surfTolerance)
303  {
304  //
305  // We are just inside or away from the
306  // surface. Accept *any* value of distance.
307  //
308  distance = s1;
309  return true;
310  }
311  }
312  }
313  else
314  distFromSurface = s1;
315 
316  //
317  // Accept positive distances
318  //
319  if (s1 > 0)
320  {
321  distance = s1;
322  return true;
323  }
324  }
325  }
326 
327  if (nside == 1) return false;
328 
329  //
330  // Well, try the second hit
331  //
332  hit = p + s2 * v;
333 
334  if (PointOnCone(hit, normSign, p, v, normal))
335  {
336  //
337  // Good intersection! What about the normal?
338  //
339  if (normSign * v.Dot(normal) > 0)
340  {
341  double pr = p.Perp();
342  if (pr < DBL_MIN) pr = DBL_MIN;
343  UVector3 pNormal(rNorm * p.x / pr, rNorm * p.y / pr, zNorm);
344  if (normSign * v.Dot(pNormal) > 0)
345  {
346  double distOutside2;
347  distFromSurface = -normSign * DistanceAway(p, false, distOutside2);
348  if (distOutside2 < surfTolerance * surfTolerance)
349  {
350  if (distFromSurface > -surfTolerance)
351  {
352  distance = s2;
353  return true;
354  }
355  }
356  }
357  else
358  distFromSurface = s2;
359 
360  if (s2 > 0)
361  {
362  distance = s2;
363  return true;
364  }
365  }
366  }
367 
368  //
369  // Better luck next time
370  //
371  return false;
372 }
373 
374 
375 double UPolyconeSide::Safety(const UVector3& p, bool outgoing)
376 {
377  double normSign = outgoing ? -1 : +1;
378  double distFrom, distOut2;
379 
380  //
381  // We have two tries for each hemisphere. Try the closest first.
382  //
383  distFrom = normSign * DistanceAway(p, false, distOut2);
384  if (distFrom > -0.5 * VUSolid::Tolerance())
385  {
386  //
387  // Good answer
388  //
389  if (distOut2 > 0)
390  return std::sqrt(distFrom * distFrom + distOut2);
391  else
392  return std::fabs(distFrom);
393  }
394 
395  //
396  // Try second side.
397  //
398  distFrom = normSign * DistanceAway(p, true, distOut2);
399  if (distFrom > -0.5 * VUSolid::Tolerance())
400  {
401 
402  if (distOut2 > 0)
403  return std::sqrt(distFrom * distFrom + distOut2);
404  else
405  return std::fabs(distFrom);
406  }
407 
408  return UUtils::kInfinity;
409 }
410 
411 
412 //
413 // Inside
414 //
416  double atolerance,
417  double* bestDistance)
418 {
419  //
420  // Check both sides
421  //
422  double distFrom[2], distOut2[2], dist2[2];
423  double edgeRZnorm[2];
424 
425  distFrom[0] = DistanceAway(p, false, distOut2[0], edgeRZnorm);
426  distFrom[1] = DistanceAway(p, true, distOut2[1], edgeRZnorm + 1);
427 
428  dist2[0] = distFrom[0] * distFrom[0] + distOut2[0];
429  dist2[1] = distFrom[1] * distFrom[1] + distOut2[1];
430 
431  //
432  // Who's closest?
433  //
434  int i = std::fabs(dist2[0]) < std::fabs(dist2[1]) ? 0 : 1;
435 
436  *bestDistance = std::sqrt(dist2[i]); // could sqrt be removed?
437 
438  //
439  // Okay then, inside or out?
440  //
441  if ((std::fabs(edgeRZnorm[i]) < atolerance)
442  && (distOut2[i] < atolerance * atolerance))
443  return VUSolid::eSurface;
444  else if (edgeRZnorm[i] < 0)
445  return VUSolid::eInside;
446  else
447  return VUSolid::eOutside;
448 }
449 
450 
451 //
452 // Normal
453 //
455  double* bestDistance)
456 {
457  if (p == UVector3(0., 0., 0.))
458  {
459  return p;
460  }
461 
462  double dFrom, dOut2;
463 
464  dFrom = DistanceAway(p, false, dOut2);
465 
466  *bestDistance = std::sqrt(dFrom * dFrom + dOut2);
467 
468  double rds = p.Perp();
469  if (rds != 0.)
470  {
471  return UVector3(rNorm * p.x / rds, rNorm * p.y / rds, zNorm);
472  }
473  return UVector3(0., 0., zNorm).Unit();
474 }
475 
476 
477 //
478 // Extent
479 //
480 double UPolyconeSide::Extent(const UVector3 axis)
481 {
482  if (axis.Perp2() < DBL_MIN)
483  {
484  //
485  // Special case
486  //
487  return axis.z < 0 ? -cone->ZLo() : cone->ZHi();
488  }
489 
490  //
491  // Is the axis pointing inside our phi gap?
492  //
493  if (phiIsOpen)
494  {
495  double phi = GetPhi(axis);
496  while (phi < startPhi) phi += 2 * UUtils::kPi;
497 
498  if (phi > deltaPhi + startPhi)
499  {
500  //
501  // Yeah, looks so. Make four three vectors defining the phi
502  // opening
503  //
504  double cosP = std::cos(startPhi), sinP = std::sin(startPhi);
505  UVector3 a(r[0]*cosP, r[0]*sinP, z[0]);
506  UVector3 b(r[1]*cosP, r[1]*sinP, z[1]);
507  cosP = std::cos(startPhi + deltaPhi);
508  sinP = std::sin(startPhi + deltaPhi);
509  UVector3 c(r[0]*cosP, r[0]*sinP, z[0]);
510  UVector3 d(r[1]*cosP, r[1]*sinP, z[1]);
511 
512  double ad = axis.Dot(a),
513  bd = axis.Dot(b),
514  cd = axis.Dot(c),
515  dd = axis.Dot(d);
516 
517  if (bd > ad) ad = bd;
518  if (cd > ad) ad = cd;
519  if (dd > ad) ad = dd;
520 
521  return ad;
522  }
523  }
524 
525  //
526  // Check either end
527  //
528  double aPerp = axis.Perp();
529 
530  double a = aPerp * r[0] + axis.z * z[0];
531  double b = aPerp * r[1] + axis.z * z[1];
532 
533  if (b > a) a = b;
534 
535  return a;
536 }
537 
538 
539 
540 //
541 // CalculateExtent
542 //
543 // See notes in UVCSGface
544 //
545 
546 /*
547 void UPolyconeSide::CalculateExtent( const EAxisType axis,
548  const UVoxelLimits &voxelLimit,
549  const UAffineTransform &transform,
550  USolidExtentList &extentList )
551 {
552  UClippablePolygon polygon;
553 
554  //
555  // Here we will approximate (ala UCons) and divide our conical section
556  // into segments, like UPolyhedra. When doing so, the radius
557  // is extented far enough such that the segments always lie
558  // just outside the surface of the conical section we are
559  // approximating.
560  //
561 
562  //
563  // Choose phi size of our segment(s) based on constants as
564  // defined in meshdefs.hh
565  //
566  int numPhi = (int)(deltaPhi/UUtils::kMeshAngleDefault) + 1;
567  if (numPhi < UUtils::kMinMeshSections)
568  numPhi = UUtils::kMinMeshSections;
569  else if (numPhi > UUtils::kMaxMeshSections)
570  numPhi = UUtils::kMaxMeshSections;
571 
572  double sigPhi = deltaPhi/numPhi;
573 
574  //
575  // Determine radius factor to keep segments outside
576  //
577  double rFudge = 1.0/std::cos(0.5*sigPhi);
578 
579  //
580  // Decide which radius to use on each end of the side,
581  // and whether a transition mesh is required
582  //
583  // {r0,z0} - Beginning of this side
584  // {r1,z1} - Ending of this side
585  // {r2,z0} - Beginning of transition piece connecting previous
586  // side (and ends at beginning of this side)
587  //
588  // So, order is 2 --> 0 --> 1.
589  // -------
590  //
591  // r2 < 0 indicates that no transition piece is required
592  //
593  double r0, r1, r2, z0, z1;
594 
595  r2 = -1; // By default: no transition piece
596 
597  if (rNorm < -DBL_MIN)
598  {
599  //
600  // This side faces *inward*, and so our mesh has
601  // the same radius
602  //
603  r1 = r[1];
604  z1 = z[1];
605  z0 = z[0];
606  r0 = r[0];
607 
608  r2 = -1;
609 
610  if (prevZS > DBL_MIN)
611  {
612  //
613  // The previous side is facing outwards
614  //
615  if ( prevRS*zS - prevZS*rS > 0 )
616  {
617  //
618  // Transition was convex: build transition piece
619  //
620  if (r[0] > DBL_MIN) r2 = r[0]*rFudge;
621  }
622  else
623  {
624  //
625  // Transition was concave: short this side
626  //
627  FindLineIntersect( z0, r0, zS, rS,
628  z0, r0*rFudge, prevZS, prevRS*rFudge, z0, r0 );
629  }
630  }
631 
632  if ( nextZS > DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
633  {
634  //
635  // The next side is facing outwards, forming a
636  // concave transition: short this side
637  //
638  FindLineIntersect( z1, r1, zS, rS,
639  z1, r1*rFudge, nextZS, nextRS*rFudge, z1, r1 );
640  }
641  }
642  else if (rNorm > DBL_MIN)
643  {
644  //
645  // This side faces *outward* and is given a boost to
646  // it radius
647  //
648  r0 = r[0]*rFudge;
649  z0 = z[0];
650  r1 = r[1]*rFudge;
651  z1 = z[1];
652 
653  if (prevZS < -DBL_MIN)
654  {
655  //
656  // The previous side is facing inwards
657  //
658  if ( prevRS*zS - prevZS*rS > 0 )
659  {
660  //
661  // Transition was convex: build transition piece
662  //
663  if (r[0] > DBL_MIN) r2 = r[0];
664  }
665  else
666  {
667  //
668  // Transition was concave: short this side
669  //
670  FindLineIntersect( z0, r0, zS, rS*rFudge,
671  z0, r[0], prevZS, prevRS, z0, r0 );
672  }
673  }
674 
675  if ( nextZS < -DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
676  {
677  //
678  // The next side is facing inwards, forming a
679  // concave transition: short this side
680  //
681  FindLineIntersect( z1, r1, zS, rS*rFudge,
682  z1, r[1], nextZS, nextRS, z1, r1 );
683  }
684  }
685  else
686  {
687  //
688  // This side is perpendicular to the z axis (is a disk)
689  //
690  // Whether or not r0 needs a rFudge factor depends
691  // on the normal of the previous edge. Similar with r1
692  // and the next edge. No transition piece is required.
693  //
694  r0 = r[0];
695  r1 = r[1];
696  z0 = z[0];
697  z1 = z[1];
698 
699  if (prevZS > DBL_MIN) r0 *= rFudge;
700  if (nextZS > DBL_MIN) r1 *= rFudge;
701  }
702 
703  //
704  // Loop
705  //
706  double phi = startPhi,
707  cosPhi = std::cos(phi),
708  sinPhi = std::sin(phi);
709 
710  UVector3 v0( r0*cosPhi, r0*sinPhi, z0 ),
711  v1( r1*cosPhi, r1*sinPhi, z1 ),
712  v2, w0, w1, w2;
713  transform.ApplyPointTransform( v0 );
714  transform.ApplyPointTransform( v1 );
715 
716  if (r2 >= 0)
717  {
718  v2 = UVector3( r2*cosPhi, r2*sinPhi, z0 );
719  transform.ApplyPointTransform( v2 );
720  }
721 
722  do
723  {
724  phi += sigPhi;
725  if (numPhi == 1) phi = startPhi+deltaPhi; // Try to avoid roundoff
726  cosPhi = std::cos(phi),
727  sinPhi = std::sin(phi);
728 
729  w0 = UVector3( r0*cosPhi, r0*sinPhi, z0 );
730  w1 = UVector3( r1*cosPhi, r1*sinPhi, z1 );
731  transform.ApplyPointTransform( w0 );
732  transform.ApplyPointTransform( w1 );
733 
734  UVector3 deltaV = r0 > r1 ? w0-v0 : w1-v1;
735 
736  //
737  // Build polygon, taking special care to keep the vertices
738  // in order
739  //
740  polygon.ClearAllVertices();
741 
742  polygon.AddVertexInOrder( v0 );
743  polygon.AddVertexInOrder( v1 );
744  polygon.AddVertexInOrder( w1 );
745  polygon.AddVertexInOrder( w0 );
746 
747  //
748  // Get extent
749  //
750  if (polygon.PartialClip( voxelLimit, axis ))
751  {
752  //
753  // Get Dot product of normal with target axis
754  //
755  polygon.SetNormal( deltaV.Cross(v1-v0).Unit() );
756 
757  extentList.AddSurface( polygon );
758  }
759 
760  if (r2 >= 0)
761  {
762  //
763  // Repeat, for transition piece
764  //
765  w2 = UVector3( r2*cosPhi, r2*sinPhi, z0 );
766  transform.ApplyPointTransform( w2 );
767 
768  polygon.ClearAllVertices();
769 
770  polygon.AddVertexInOrder( v2 );
771  polygon.AddVertexInOrder( v0 );
772  polygon.AddVertexInOrder( w0 );
773  polygon.AddVertexInOrder( w2 );
774 
775  if (polygon.PartialClip( voxelLimit, axis ))
776  {
777  polygon.SetNormal( deltaV.Cross(v0-v2).Unit() );
778 
779  extentList.AddSurface( polygon );
780  }
781 
782  v2 = w2;
783  }
784 
785  //
786  // Next vertex
787  //
788  v0 = w0;
789  v1 = w1;
790  } while( --numPhi > 0 );
791 
792  //
793  // We are almost done. But, it is important that we leave no
794  // gaps in the surface of our solid. By using rFudge, however,
795  // we've done exactly that, if we have a phi segment.
796  // Add two additional faces if necessary
797  //
798  if (phiIsOpen && rNorm > DBL_MIN)
799  {
800  cosPhi = std::cos(startPhi);
801  sinPhi = std::sin(startPhi);
802 
803  UVector3 a0( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
804  a1( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
805  b0( r0*cosPhi, r0*sinPhi, z[0] ),
806  b1( r1*cosPhi, r1*sinPhi, z[1] );
807 
808  transform.ApplyPointTransform( a0 );
809  transform.ApplyPointTransform( a1 );
810  transform.ApplyPointTransform( b0 );
811  transform.ApplyPointTransform( b1 );
812 
813  polygon.ClearAllVertices();
814 
815  polygon.AddVertexInOrder( a0 );
816  polygon.AddVertexInOrder( a1 );
817  polygon.AddVertexInOrder( b0 );
818  polygon.AddVertexInOrder( b1 );
819 
820  if (polygon.PartialClip( voxelLimit , axis))
821  {
822  UVector3 normal( sinPhi, -cosPhi, 0 );
823  polygon.SetNormal( transform.TransformAxis( normal ) );
824 
825  extentList.AddSurface( polygon );
826  }
827 
828  cosPhi = std::cos(startPhi+deltaPhi);
829  sinPhi = std::sin(startPhi+deltaPhi);
830 
831  a0 = UVector3( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
832  a1 = UVector3( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
833  b0 = UVector3( r0*cosPhi, r0*sinPhi, z[0] ),
834  b1 = UVector3( r1*cosPhi, r1*sinPhi, z[1] );
835  transform.ApplyPointTransform( a0 );
836  transform.ApplyPointTransform( a1 );
837  transform.ApplyPointTransform( b0 );
838  transform.ApplyPointTransform( b1 );
839 
840  polygon.ClearAllVertices();
841 
842  polygon.AddVertexInOrder( a0 );
843  polygon.AddVertexInOrder( a1 );
844  polygon.AddVertexInOrder( b0 );
845  polygon.AddVertexInOrder( b1 );
846 
847  if (polygon.PartialClip( voxelLimit, axis ))
848  {
849  UVector3 normal( -sinPhi, cosPhi, 0 );
850  polygon.SetNormal( transform.TransformAxis( normal ) );
851 
852  extentList.AddSurface( polygon );
853  }
854  }
855 
856  return;
857 }
858 */
859 
860 //
861 // GetPhi
862 //
863 // Calculate Phi for a given 3-vector (point), if not already cached for the
864 // same point, in the attempt to avoid consecutive computation of the same
865 // quantity
866 //
868 {
869  double val = 0.;
870 
871  val = p.Phi();
872 
873  return val;
874 }
875 
876 //
877 // DistanceAway
878 //
879 // Calculate distance of a point from our conical surface, including the effect
880 // of any phi segmentation
881 //
882 // Arguments:
883 // p - (in) Point to check
884 // opposite - (in) If true, check opposite hemisphere (see below)
885 // distOutside - (out) Additional distance outside the edges of the surface
886 // edgeRZnorm - (out) if negative, point is inside
887 //
888 // return value = distance from the conical plane, if extrapolated beyond edges,
889 // signed by whether the point is in inside or outside the shape
890 //
891 // Notes:
892 // * There are two answers, depending on which hemisphere is considered.
893 //
895  bool opposite,
896  double& distOutside2,
897  double* edgeRZnorm)
898 {
899  //
900  // Convert our point to r and z
901  //
902  double rx = p.Perp(), zx = p.z;
903 
904  //
905  // Change sign of r if opposite says we should
906  //
907  if (opposite) rx = -rx;
908 
909  //
910  // Calculate return value
911  //
912  double deltaR = rx - r[0], deltaZ = zx - z[0];
913  double answer = deltaR * rNorm + deltaZ * zNorm;
914 
915  //
916  // Are we off the surface in r,z space?
917  //
918  double q = deltaR * rS + deltaZ * zS;
919  if (q < 0)
920  {
921  distOutside2 = q * q;
922  if (edgeRZnorm) *edgeRZnorm = deltaR * rNormEdge[0] + deltaZ * zNormEdge[0];
923  }
924  else if (q > length)
925  {
926  distOutside2 = UUtils::sqr(q - length);
927  if (edgeRZnorm)
928  {
929  deltaR = rx - r[1];
930  deltaZ = zx - z[1];
931  *edgeRZnorm = deltaR * rNormEdge[1] + deltaZ * zNormEdge[1];
932  }
933  }
934  else
935  {
936  distOutside2 = 0;
937  if (edgeRZnorm) *edgeRZnorm = answer;
938  }
939 
940  if (phiIsOpen)
941  {
942  //
943  // Finally, check phi
944  //
945  double phi = GetPhi(p);
946  while (phi < startPhi) phi += 2 * UUtils::kPi;
947 
948  if (phi > startPhi + deltaPhi)
949  {
950  //
951  // Oops. Are we closer to the start phi or end phi?
952  //
953  double d1 = phi - startPhi - deltaPhi;
954  while (phi > startPhi) phi -= 2 * UUtils::kPi;
955  double d2 = startPhi - phi;
956 
957  if (d2 < d1) d1 = d2;
958 
959  //
960  // Add result to our distance
961  //
962  double dist = d1 * rx;
963 
964  distOutside2 += dist * dist;
965  if (edgeRZnorm)
966  {
967  *edgeRZnorm = std::max(std::fabs(*edgeRZnorm), std::fabs(dist));
968  }
969  }
970  }
971 
972  return answer;
973 }
974 
975 
976 //
977 // PointOnCone
978 //
979 // Decide if a point is on a cone and return normal if it is
980 //
982  double normSign,
983  const UVector3& p,
984  const UVector3& v,
985  UVector3& normal)
986 {
987  double rx = hit.Perp();
988  //
989  // Check radial/z extent, as appropriate
990  //
991  if (!cone->HitOn(rx, hit.z)) return false;
992 
993  if (phiIsOpen)
994  {
995  double phiTolerant = 2.0 * VUSolid::Tolerance() / (rx + VUSolid::Tolerance());
996  //
997  // Check phi segment. Here we have to be careful
998  // to use the standard method consistent with
999  // PolyPhiFace. See PolyPhiFace::InsideEdgesExact
1000  //
1001  double phi = GetPhi(hit);
1002  while (phi < startPhi - phiTolerant) phi += 2 * UUtils::kPi;
1003 
1004  if (phi > startPhi + deltaPhi + phiTolerant) return false;
1005 
1006  if (phi > startPhi + deltaPhi - phiTolerant)
1007  {
1008  //
1009  // Exact treatment
1010  //
1011  UVector3 qx = p + v;
1012  UVector3 qa = qx - corners[2],
1013  qb = qx - corners[3];
1014  UVector3 qacb = qa.Cross(qb);
1015 
1016  if (normSign * qacb.Dot(v) < 0) return false;
1017  }
1018  else if (phi < phiTolerant)
1019  {
1020  UVector3 qx = p + v;
1021  UVector3 qa = qx - corners[1],
1022  qb = qx - corners[0];
1023  UVector3 qacb = qa.Cross(qb);
1024 
1025  if (normSign * qacb.Dot(v) < 0) return false;
1026  }
1027  }
1028 
1029  //
1030  // We have a good hit! Calculate normal
1031  //
1032  if (rx < DBL_MIN)
1033  normal = UVector3(0, 0, zNorm < 0 ? -1 : 1);
1034  else
1035  normal = UVector3(rNorm * hit.x / rx, rNorm * hit.y / rx, zNorm);
1036  return true;
1037 }
1038 
1039 
1040 //
1041 // FindLineIntersect
1042 //
1043 // Decide the point at which two 2-dimensional lines intersect
1044 //
1045 // Equation of line: x = x1 + s*tx1
1046 // y = y1 + s*ty1
1047 //
1048 // It is assumed that the lines are *not* parallel
1049 //
1050 void UPolyconeSide::FindLineIntersect(double x1, double y1,
1051  double tx1, double ty1,
1052  double x2, double y2,
1053  double tx2, double ty2,
1054  double& x, double& y)
1055 {
1056  //
1057  // The solution is a simple linear equation
1058  //
1059  double deter = tx1 * ty2 - tx2 * ty1;
1060 
1061  double s1 = ((x2 - x1) * ty2 - tx2 * (y2 - y1)) / deter;
1062  double s2 = ((x2 - x1) * ty1 - tx1 * (y2 - y1)) / deter;
1063 
1064  //
1065  // We want the answer to not depend on which order the
1066  // lines were specified. Take average.
1067  //
1068  x = 0.5 * (x1 + s1 * tx1 + x2 + s2 * tx2);
1069  y = 0.5 * (y1 + s1 * ty1 + y2 + s2 * ty2);
1070 }
1071 
1072 //
1073 // Calculate surface area for GetPointOnSurface()
1074 //
1076 {
1077  if (fSurfaceArea == 0)
1078  {
1079  fSurfaceArea = (r[0] + r[1]) * std::sqrt(UUtils::sqr(r[0] - r[1]) + UUtils::sqr(z[0] - z[1]));
1080  fSurfaceArea *= 0.5 * (deltaPhi);
1081  }
1082  return fSurfaceArea;
1083 }
1084 
1085 //
1086 // GetPointOnFace
1087 //
1089 {
1090  double x, y, zz;
1091  double rr, phi, dz, dr;
1092  dr = r[1] - r[0];
1093  dz = z[1] - z[0];
1094  phi = startPhi + deltaPhi * UUtils::Random();
1095  rr = r[0] + dr * UUtils::Random();
1096 
1097  x = rr * std::cos(phi);
1098  y = rr * std::sin(phi);
1099 
1100  // PolyconeSide has a Ring Form
1101  //
1102  if (dz == 0.)
1103  {
1104  zz = z[0];
1105  }
1106  else
1107  {
1108  if (dr == 0.) // PolyconeSide has a Tube Form
1109  {
1110  zz = z[0] + dz * UUtils::Random();
1111  }
1112  else
1113  {
1114  zz = z[0] + (rr - r[0]) * dz / dr;
1115  }
1116  }
1117 
1118  return UVector3(x, y, zz);
1119 }
static const G4double d1
UVector3 GetPointOnFace()
void CopyStuff(const UPolyconeSide &source)
double Phi() const
Definition: UVector3.cc:64
double DistanceAway(const UVector3 &p, bool opposite, double &distOutside2, double *rzNorm=0)
UVector3 Cross(const UVector3 &) const
Definition: UVector3.hh:262
double Extent(const UVector3 axis)
double zNormEdge[2]
UPolyconeSide & operator=(const UPolyconeSide &source)
VUSolid::EnumInside Inside(const UVector3 &p, double tolerance, double *bestDistance)
static const G4double cd
G4double a
Definition: TRTMaterials.hh:39
double GetPhi(const UVector3 &p)
static double Tolerance()
Definition: VUSolid.hh:127
double x
Definition: UVector3.hh:136
bool HitOn(const double r, const double z)
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
static const double kInfinity
Definition: UUtils.hh:53
UVector3 Normal(const UVector3 &p, double *bestDistance)
double SurfaceArea()
EnumInside
Definition: VUSolid.hh:23
double Dot(const UVector3 &) const
Definition: UVector3.hh:257
T sqr(const T &x)
Definition: UUtils.hh:137
double Perp2() const
Definition: UVector3.hh:272
T max(const T t1, const T t2)
brief Return the largest of the two arguments
UIntersectingCone * cone
static const double kPi
Definition: UUtils.hh:48
double ZHi() const
UVector3 Unit() const
Definition: UVector3.cc:80
double Safety(const UVector3 &p, bool outgoing)
#define DBL_MIN
Definition: templates.hh:75
virtual ~UPolyconeSide()
UPolyconeSide(const UPolyconeSideRZ *prevRZ, const UPolyconeSideRZ *tail, const UPolyconeSideRZ *head, const UPolyconeSideRZ *nextRZ, double phiStart, double deltaPhi, bool phiIsOpen, bool isAllBehind=false)
double Perp() const
Definition: UVector3.cc:56
UVector3 * corners
double z
Definition: UVector3.hh:138
bool Distance(const UVector3 &p, const UVector3 &v, bool outgoing, double surfTolerance, double &distance, double &distFromSurface, UVector3 &normal, bool &isAllBehind)
static void FindLineIntersect(double x1, double y1, double tx1, double ty1, double x2, double y2, double tx2, double ty2, double &x, double &y)
static const G4double d2
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double ZLo() const
double y
Definition: UVector3.hh:137
int LineHitsCone(const UVector3 &p, const UVector3 &v, double &s1, double &s2)
double rNormEdge[2]
bool PointOnCone(const UVector3 &hit, double normSign, const UVector3 &p, const UVector3 &v, UVector3 &normal)