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