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