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