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