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