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