Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EllipticalTube.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: G4EllipticalTube.cc 101118 2016-11-07 09:10:59Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4EllipticalTube.cc
35 //
36 // Implementation of a CSG volume representing a tube with elliptical cross
37 // section (geant3 solid 'ELTU')
38 //
39 // --------------------------------------------------------------------
40 
41 #include "G4EllipticalTube.hh"
42 
43 #include "G4ClippablePolygon.hh"
44 #include "G4AffineTransform.hh"
45 #include "G4SolidExtentList.hh"
46 #include "G4VoxelLimits.hh"
47 #include "G4BoundingEnvelope.hh"
48 #include "meshdefs.hh"
49 
50 #include "Randomize.hh"
51 
52 #include "G4VGraphicsScene.hh"
53 #include "G4VisExtent.hh"
54 
55 #include "G4AutoLock.hh"
56 
57 namespace
58 {
59  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
60 }
61 
62 using namespace CLHEP;
63 
64 //
65 // Constructor
66 //
68  G4double theDx,
69  G4double theDy,
70  G4double theDz )
71  : G4VSolid( name ), fCubicVolume(0.), fSurfaceArea(0.),
72  fRebuildPolyhedron(false), fpPolyhedron(0)
73 {
74  halfTol = 0.5*kCarTolerance;
75 
76  dx = theDx;
77  dy = theDy;
78  dz = theDz;
79 }
80 
81 
82 //
83 // Fake default constructor - sets only member data and allocates memory
84 // for usage restricted to object persistency.
85 //
87  : G4VSolid(a), dx(0.), dy(0.), dz(0.), halfTol(0.),
88  fCubicVolume(0.), fSurfaceArea(0.),
89  fRebuildPolyhedron(false), fpPolyhedron(0)
90 {
91 }
92 
93 
94 //
95 // Destructor
96 //
98 {
99  delete fpPolyhedron; fpPolyhedron = 0;
100 }
101 
102 
103 //
104 // Copy constructor
105 //
107  : G4VSolid(rhs), dx(rhs.dx), dy(rhs.dy), dz(rhs.dz), halfTol(rhs.halfTol),
108  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
109  fRebuildPolyhedron(false), fpPolyhedron(0)
110 {
111 }
112 
113 
114 //
115 // Assignment operator
116 //
118 {
119  // Check assignment to self
120  //
121  if (this == &rhs) { return *this; }
122 
123  // Copy base class data
124  //
125  G4VSolid::operator=(rhs);
126 
127  // Copy data
128  //
129  dx = rhs.dx; dy = rhs.dy; dz = rhs.dz;
130  halfTol = rhs.halfTol;
131  fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
132  fRebuildPolyhedron = false;
133  delete fpPolyhedron; fpPolyhedron = 0;
134 
135  return *this;
136 }
137 
139 //
140 // Get bounding box
141 
143  G4ThreeVector& pMax ) const
144 {
145  pMin.set(-dx,-dy,-dz);
146  pMax.set( dx, dy, dz);
147 }
148 
150 //
151 // Calculate extent under transform and specified limit
152 
153 G4bool
155  const G4VoxelLimits& pVoxelLimit,
156  const G4AffineTransform& pTransform,
157  G4double& pMin, G4double& pMax ) const
158 {
159  G4ThreeVector bmin, bmax;
160  G4bool exist;
161 
162  // Check bounding box (bbox)
163  //
164  Extent(bmin,bmax);
165  G4BoundingEnvelope bbox(bmin,bmax);
166 #ifdef G4BBOX_EXTENT
167  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
168 #endif
169  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
170  {
171  return exist = (pMin < pMax) ? true : false;
172  }
173 
174  // Set bounding envelope (benv) and calculate extent
175  //
176  const G4int NSTEPS = 48; // number of steps for whole circle
177  G4double ang = twopi/NSTEPS;
178 
179  G4double sinHalf = std::sin(0.5*ang);
180  G4double cosHalf = std::cos(0.5*ang);
181  G4double sinStep = 2.*sinHalf*cosHalf;
182  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
183  G4double sx = dx/cosHalf;
184  G4double sy = dy/cosHalf;
185 
186  G4double sinCur = sinHalf;
187  G4double cosCur = cosHalf;
188  G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
189  for (G4int k=0; k<NSTEPS; ++k)
190  {
191  baseA[k].set(sx*cosCur,sy*sinCur,-dz);
192  baseB[k].set(sx*cosCur,sy*sinCur, dz);
193 
194  G4double sinTmp = sinCur;
195  sinCur = sinCur*cosStep + cosCur*sinStep;
196  cosCur = cosCur*cosStep - sinTmp*sinStep;
197  }
198 
199  std::vector<const G4ThreeVectorList *> polygons(2);
200  polygons[0] = &baseA;
201  polygons[1] = &baseB;
202  G4BoundingEnvelope benv(bmin,bmax,polygons);
203  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
204  return exist;
205 }
206 
207 //
208 // Inside
209 //
210 // Note that for this solid, we've decided to define the tolerant
211 // surface as that which is bounded by ellipses with axes
212 // at +/- 0.5*kCarTolerance.
213 //
215 {
216  //
217  // Check z extents: are we outside?
218  //
219  G4double absZ = std::fabs(p.z());
220  if (absZ > dz+halfTol) return kOutside;
221 
222  //
223  // Check x,y: are we outside?
224  //
225  // G4double x = p.x(), y = p.y();
226 
227  if (CheckXY(p.x(), p.y(), +halfTol) > 1.0) return kOutside;
228 
229  //
230  // We are either inside or on the surface: recheck z extents
231  //
232  if (absZ > dz-halfTol) return kSurface;
233 
234  //
235  // Recheck x,y
236  //
237  if (CheckXY(p.x(), p.y(), -halfTol) > 1.0) return kSurface;
238 
239  return kInside;
240 }
241 
242 
243 //
244 // SurfaceNormal
245 //
247 {
248  //
249  // SurfaceNormal for the point On the Surface, sum the normals on the Corners
250  //
251 
252  G4int noSurfaces=0;
253  G4ThreeVector norm, sumnorm(0.,0.,0.);
254 
255  G4double distZ = std::fabs(std::fabs(p.z()) - dz);
256 
257  G4double distR1 = CheckXY( p.x(), p.y(),+ halfTol );
258  G4double distR2 = CheckXY( p.x(), p.y(),- halfTol );
259 
260  if ( (distZ < halfTol ) && ( distR1 <= 1 ) )
261  {
262  noSurfaces++;
263  sumnorm=G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
264  }
265  if( (distR1 <= 1 ) && ( distR2 >= 1 ) )
266  {
267  noSurfaces++;
268  norm= G4ThreeVector( p.x()*dy*dy, p.y()*dx*dx, 0.0 ).unit();
269  sumnorm+=norm;
270  }
271  if ( noSurfaces == 0 )
272  {
273 #ifdef G4SPECSDEBUG
274  G4Exception("G4EllipticalTube::SurfaceNormal(p)", "GeomSolids1002",
275  JustWarning, "Point p is not on surface !?" );
276 #endif
277  norm = ApproxSurfaceNormal(p);
278  }
279  else if ( noSurfaces == 1 ) { norm = sumnorm; }
280  else { norm = sumnorm.unit(); }
281 
282  return norm;
283 }
284 
285 
286 //
287 // ApproxSurfaceNormal
288 //
290 G4EllipticalTube::ApproxSurfaceNormal( const G4ThreeVector& p ) const
291 {
292  //
293  // Which of the three surfaces are we closest to (approximatively)?
294  //
295  G4double distZ = std::fabs(p.z()) - dz;
296 
297  G4double rxy = CheckXY( p.x(), p.y() );
298  G4double distR2 = (rxy < DBL_MIN) ? DBL_MAX : 1.0/rxy;
299 
300  //
301  // Closer to z?
302  //
303  if (distZ*distZ < distR2)
304  {
305  return G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
306  }
307 
308  //
309  // Closer to x/y
310  //
311  return G4ThreeVector( p.x()*dy*dy, p.y()*dx*dx, 0.0 ).unit();
312 }
313 
314 
315 //
316 // DistanceToIn(p,v)
317 //
318 // Unlike DistanceToOut(p,v), it is possible for the trajectory
319 // to miss. The geometric calculations here are quite simple.
320 // More difficult is the logic required to prevent particles
321 // from sneaking (or leaking) between the elliptical and end
322 // surfaces.
323 //
324 // Keep in mind that the true distance is allowed to be
325 // negative if the point is currently on the surface. For oblique
326 // angles, it can be very negative.
327 //
329  const G4ThreeVector& v ) const
330 {
331  //
332  // Check z = -dz planer surface
333  //
334  G4double sigz = p.z()+dz;
335 
336  if (sigz < halfTol)
337  {
338  //
339  // We are "behind" the shape in z, and so can
340  // potentially hit the rear face. Correct direction?
341  //
342  if (v.z() <= 0)
343  {
344  //
345  // As long as we are far enough away, we know we
346  // can't intersect
347  //
348  if (sigz < 0) return kInfinity;
349 
350  //
351  // Otherwise, we don't intersect unless we are
352  // on the surface of the ellipse
353  //
354  if (CheckXY(p.x(),p.y(),-halfTol) <= 1.0) return kInfinity;
355  }
356  else
357  {
358  //
359  // How far?
360  //
361  G4double q = -sigz/v.z();
362 
363  //
364  // Where does that place us?
365  //
366  G4double xi = p.x() + q*v.x(),
367  yi = p.y() + q*v.y();
368 
369  //
370  // Is this on the surface (within ellipse)?
371  //
372  if (CheckXY(xi,yi) <= 1.0)
373  {
374  //
375  // Yup. Return q, unless we are on the surface
376  //
377  return (sigz < -halfTol) ? q : 0;
378  }
379  else if (xi*dy*dy*v.x() + yi*dx*dx*v.y() >= 0)
380  {
381  //
382  // Else, if we are traveling outwards, we know
383  // we must miss
384  //
385  return kInfinity;
386  }
387  }
388  }
389 
390  //
391  // Check z = +dz planer surface
392  //
393  sigz = p.z() - dz;
394 
395  if (sigz > -halfTol)
396  {
397  if (v.z() >= 0)
398  {
399  if (sigz > 0) return kInfinity;
400  if (CheckXY(p.x(),p.y(),-halfTol) <= 1.0) return kInfinity;
401  }
402  else {
403  G4double q = -sigz/v.z();
404 
405  G4double xi = p.x() + q*v.x(),
406  yi = p.y() + q*v.y();
407 
408  if (CheckXY(xi,yi) <= 1.0)
409  {
410  return (sigz > -halfTol) ? q : 0;
411  }
412  else if (xi*dy*dy*v.x() + yi*dx*dx*v.y() >= 0)
413  {
414  return kInfinity;
415  }
416  }
417  }
418 
419  //
420  // Check intersection with the elliptical tube
421  //
422  G4double q[2];
423  G4int n = IntersectXY( p, v, q );
424 
425  if (n==0) return kInfinity;
426 
427  //
428  // Is the original point on the surface?
429  //
430  if (std::fabs(p.z()) < dz+halfTol) {
431  if (CheckXY( p.x(), p.y(), halfTol ) < 1.0)
432  {
433  //
434  // Well, yes, but are we traveling inwards at this point?
435  //
436  if (p.x()*dy*dy*v.x() + p.y()*dx*dx*v.y() < 0) return 0;
437  }
438  }
439 
440  //
441  // We are now certain that point p is not on the surface of
442  // the solid (and thus std::fabs(q[0]) > halfTol).
443  // Return kInfinity if the intersection is "behind" the point.
444  //
445  if (q[0] < 0) return kInfinity;
446 
447  //
448  // Check to see if we intersect the tube within
449  // dz, but only when we know it might miss
450  //
451  G4double zi = p.z() + q[0]*v.z();
452 
453  if (v.z() < 0)
454  {
455  if (zi < -dz) return kInfinity;
456  }
457  else if (v.z() > 0)
458  {
459  if (zi > +dz) return kInfinity;
460  }
461 
462  return q[0];
463 }
464 
465 
466 //
467 // DistanceToIn(p)
468 //
469 // The distance from a point to an ellipse (in 2 dimensions) is a
470 // surprisingly complicated quadric expression (this is easy to
471 // appreciate once one understands that there may be up to
472 // four lines normal to the ellipse intersecting any point). To
473 // solve it exactly would be rather time consuming. This method,
474 // however, is supposed to be a quick check, and is allowed to be an
475 // underestimate.
476 //
477 // So, I will use the following underestimate of the distance
478 // from an outside point to an ellipse. First: find the intersection "A"
479 // of the line from the origin to the point with the ellipse.
480 // Find the line passing through "A" and tangent to the ellipse
481 // at A. The distance of the point p from the ellipse will be approximated
482 // as the distance to this line.
483 //
485 {
486  if (CheckXY( p.x(), p.y(), +halfTol ) < 1.0)
487  {
488  //
489  // We are inside or on the surface of the
490  // elliptical cross section in x/y. Check z
491  //
492  if (p.z() < -dz-halfTol)
493  return -p.z()-dz;
494  else if (p.z() > dz+halfTol)
495  return p.z()-dz;
496  else
497  return 0; // On any surface here (or inside)
498  }
499 
500  //
501  // Find point on ellipse
502  //
503  G4double qnorm = CheckXY( p.x(), p.y() );
504  if (qnorm < DBL_MIN) return 0; // This should never happen
505 
506  G4double q = 1.0/std::sqrt(qnorm);
507 
508  G4double xe = q*p.x(), ye = q*p.y();
509 
510  //
511  // Get tangent to ellipse
512  //
513  G4double tx = -ye*dx*dx, ty = +xe*dy*dy;
514  G4double tnorm = std::sqrt( tx*tx + ty*ty );
515 
516  //
517  // Calculate distance
518  //
519  G4double distR = ( (p.x()-xe)*ty - (p.y()-ye)*tx )/tnorm;
520 
521  //
522  // Add the result in quadrature if we are, in addition,
523  // outside the z bounds of the shape
524  //
525  // We could save some time by returning the maximum rather
526  // than the quadrature sum
527  //
528  if (p.z() < -dz)
529  return std::sqrt( (p.z()+dz)*(p.z()+dz) + distR*distR );
530  else if (p.z() > dz)
531  return std::sqrt( (p.z()-dz)*(p.z()-dz) + distR*distR );
532 
533  return distR;
534 }
535 
536 
537 //
538 // DistanceToOut(p,v)
539 //
540 // This method can be somewhat complicated for a general shape.
541 // For a convex one, like this, there are several simplifications,
542 // the most important of which is that one can treat the surfaces
543 // as infinite in extent when deciding if the p is on the surface.
544 //
546  const G4ThreeVector& v,
547  const G4bool calcNorm,
548  G4bool *validNorm,
549  G4ThreeVector *norm ) const
550 {
551  //
552  // Our normal is always valid
553  //
554  if (calcNorm) { *validNorm = true; }
555 
556  G4double sBest = kInfinity;
557  G4ThreeVector nBest(0,0,0);
558 
559  //
560  // Might we intersect the -dz surface?
561  //
562  if (v.z() < 0)
563  {
564  static const G4ThreeVector normHere(0.0,0.0,-1.0);
565  //
566  // Yup. What distance?
567  //
568  sBest = -(p.z()+dz)/v.z();
569 
570  //
571  // Are we on the surface? If so, return zero
572  //
573  if (p.z() < -dz+halfTol)
574  {
575  if (calcNorm) { *norm = normHere; }
576  return 0;
577  }
578  else
579  {
580  nBest = normHere;
581  }
582  }
583 
584  //
585  // How about the +dz surface?
586  //
587  if (v.z() > 0)
588  {
589  static const G4ThreeVector normHere(0.0,0.0,+1.0);
590  //
591  // Yup. What distance?
592  //
593  G4double q = (dz-p.z())/v.z();
594 
595  //
596  // Are we on the surface? If so, return zero
597  //
598  if (p.z() > +dz-halfTol)
599  {
600  if (calcNorm) { *norm = normHere; }
601  return 0;
602  }
603 
604  //
605  // Best so far?
606  //
607  if (q < sBest) { sBest = q; nBest = normHere; }
608  }
609 
610  //
611  // Check furthest intersection with ellipse
612  //
613  G4double q[2];
614  G4int n = IntersectXY( p, v, q );
615 
616  if (n == 0)
617  {
618  if (sBest == kInfinity)
619  {
620  DumpInfo();
621  std::ostringstream message;
622  G4int oldprc = message.precision(16) ;
623  message << "Point p is outside !?" << G4endl
624  << "Position:" << G4endl
625  << " p.x() = " << p.x()/mm << " mm" << G4endl
626  << " p.y() = " << p.y()/mm << " mm" << G4endl
627  << " p.z() = " << p.z()/mm << " mm" << G4endl
628  << "Direction:" << G4endl << G4endl
629  << " v.x() = " << v.x() << G4endl
630  << " v.y() = " << v.y() << G4endl
631  << " v.z() = " << v.z() << G4endl
632  << "Proposed distance :" << G4endl
633  << " snxt = " << sBest/mm << " mm";
634  message.precision(oldprc) ;
635  G4Exception( "G4EllipticalTube::DistanceToOut(p,v,...)",
636  "GeomSolids1002", JustWarning, message);
637  }
638  if (calcNorm) { *norm = nBest; }
639  return sBest;
640  }
641  else if (q[n-1] > sBest)
642  {
643  if (calcNorm) { *norm = nBest; }
644  return sBest;
645  }
646  sBest = q[n-1];
647 
648  //
649  // Intersection with ellipse. Get normal at intersection point.
650  //
651  if (calcNorm)
652  {
653  G4ThreeVector ip = p + sBest*v;
654  *norm = G4ThreeVector( ip.x()*dy*dy, ip.y()*dx*dx, 0.0 ).unit();
655  }
656 
657  //
658  // Do we start on the surface?
659  //
660  if (CheckXY( p.x(), p.y(), -halfTol ) > 1.0)
661  {
662  //
663  // Well, yes, but are we traveling outwards at this point?
664  //
665  if (p.x()*dy*dy*v.x() + p.y()*dx*dx*v.y() > 0) return 0;
666  }
667 
668  return sBest;
669 }
670 
671 
672 //
673 // DistanceToOut(p)
674 //
675 // See DistanceToIn(p) for notes on the distance from a point
676 // to an ellipse in two dimensions.
677 //
678 // The approximation used here for a point inside the ellipse
679 // is to find the intersection with the ellipse of the lines
680 // through the point and parallel to the x and y axes. The
681 // distance of the point from the line connecting the two
682 // intersecting points is then used.
683 //
685 {
686  //
687  // We need to calculate the distances to all surfaces,
688  // and then return the smallest
689  //
690  // Check -dz and +dz surface
691  //
692  G4double sBest = dz - std::fabs(p.z());
693  if (sBest < halfTol) return 0;
694 
695  //
696  // Check elliptical surface: find intersection of
697  // line through p and parallel to x axis
698  //
699  G4double radical = 1.0 - p.y()*p.y()/dy/dy;
700  if (radical < +DBL_MIN) return 0;
701 
702  G4double xi = dx*std::sqrt( radical );
703  if (p.x() < 0) xi = -xi;
704 
705  //
706  // Do the same with y axis
707  //
708  radical = 1.0 - p.x()*p.x()/dx/dx;
709  if (radical < +DBL_MIN) return 0;
710 
711  G4double yi = dy*std::sqrt( radical );
712  if (p.y() < 0) yi = -yi;
713 
714  //
715  // Get distance from p to the line connecting
716  // these two points
717  //
718  G4double xdi = p.x() - xi,
719  ydi = yi - p.y();
720 
721  G4double normi = std::sqrt( xdi*xdi + ydi*ydi );
722  if (normi < halfTol) return 0;
723  xdi /= normi;
724  ydi /= normi;
725 
726  G4double q = 0.5*(xdi*(p.y()-yi) - ydi*(p.x()-xi));
727  if (xi*yi < 0) q = -q;
728 
729  if (q < sBest) sBest = q;
730 
731  //
732  // Return best answer
733  //
734  return sBest < halfTol ? 0 : sBest;
735 }
736 
737 
738 //
739 // IntersectXY
740 //
741 // Decide if and where the x/y trajectory hits the elliptical cross
742 // section.
743 //
744 // Arguments:
745 // p - (in) Point on trajectory
746 // v - (in) Vector along trajectory
747 // q - (out) Up to two points of intersection, where the
748 // intersection point is p + q*v, and if there are
749 // two intersections, q[0] < q[1]. May be negative.
750 // Returns:
751 // The number of intersections. If 0, the trajectory misses. If 1, the
752 // trajectory just grazes the surface.
753 //
754 // Solution:
755 // One needs to solve: ((p.x + q*v.x)/dx)**2 + ((p.y + q*v.y)/dy)**2 = 1
756 //
757 // The solution is quadratic: a*q**2 + b*q + c = 0
758 //
759 // a = (v.x/dx)**2 + (v.y/dy)**2
760 // b = 2*p.x*v.x/dx**2 + 2*p.y*v.y/dy**2
761 // c = (p.x/dx)**2 + (p.y/dy)**2 - 1
762 //
764  const G4ThreeVector &v,
765  G4double ss[2] ) const
766 {
767  G4double px = p.x(), py = p.y();
768  G4double vx = v.x(), vy = v.y();
769 
770  G4double a = (vx/dx)*(vx/dx) + (vy/dy)*(vy/dy);
771  G4double b = 2.0*( px*vx/dx/dx + py*vy/dy/dy );
772  G4double c = (px/dx)*(px/dx) + (py/dy)*(py/dy) - 1.0;
773 
774  if (a < DBL_MIN) return 0; // Trajectory parallel to z axis
775 
776  G4double radical = b*b - 4*a*c;
777 
778  if (radical < -DBL_MIN) return 0; // No solution
779 
780  if (radical < DBL_MIN)
781  {
782  //
783  // Grazes surface
784  //
785  ss[0] = -b/a/2.0;
786  return 1;
787  }
788 
789  radical = std::sqrt(radical);
790 
791  G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
792  G4double sa = q/a;
793  G4double sb = c/q;
794  if (sa < sb) { ss[0] = sa; ss[1] = sb; } else { ss[0] = sb; ss[1] = sa; }
795  return 2;
796 }
797 
798 
799 //
800 // GetEntityType
801 //
803 {
804  return G4String("G4EllipticalTube");
805 }
806 
807 
808 //
809 // Make a clone of the object
810 //
812 {
813  return new G4EllipticalTube(*this);
814 }
815 
816 
817 //
818 // GetCubicVolume
819 //
821 {
822  if(fCubicVolume != 0.) {;}
823  else { fCubicVolume = G4VSolid::GetCubicVolume(); }
824  return fCubicVolume;
825 }
826 
827 //
828 // GetSurfaceArea
829 //
831 {
832  if(fSurfaceArea != 0.) {;}
833  else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
834  return fSurfaceArea;
835 }
836 
837 //
838 // Stream object contents to an output stream
839 //
840 std::ostream& G4EllipticalTube::StreamInfo(std::ostream& os) const
841 {
842  G4int oldprc = os.precision(16);
843  os << "-----------------------------------------------------------\n"
844  << " *** Dump for solid - " << GetName() << " ***\n"
845  << " ===================================================\n"
846  << " Solid type: G4EllipticalTube\n"
847  << " Parameters: \n"
848  << " length Z: " << dz/mm << " mm \n"
849  << " surface equation in X and Y: \n"
850  << " (X / " << dx << ")^2 + (Y / " << dy << ")^2 = 1 \n"
851  << "-----------------------------------------------------------\n";
852  os.precision(oldprc);
853 
854  return os;
855 }
856 
857 
858 //
859 // GetPointOnSurface
860 //
861 // Randomly generates a point on the surface,
862 // with ~ uniform distribution across surface.
863 //
865 {
866  G4double xRand, yRand, zRand, phi, cosphi, sinphi, zArea, cArea,p, chose;
867 
868  phi = G4RandFlat::shoot(0., 2.*pi);
869  cosphi = std::cos(phi);
870  sinphi = std::sin(phi);
871 
872  // the ellipse perimeter from: "http://mathworld.wolfram.com/Ellipse.html"
873  // m = (dx - dy)/(dx + dy);
874  // k = 1.+1./4.*m*m+1./64.*sqr(m)*sqr(m)+1./256.*sqr(m)*sqr(m)*sqr(m);
875  // p = pi*(a+b)*k;
876 
877  // perimeter below from "http://www.efunda.com/math/areas/EllipseGen.cfm"
878 
879  p = 2.*pi*std::sqrt(0.5*(dx*dx+dy*dy));
880 
881  cArea = 2.*dz*p;
882  zArea = pi*dx*dy;
883 
884  xRand = dx*cosphi;
885  yRand = dy*sinphi;
886  zRand = G4RandFlat::shoot(dz, -1.*dz);
887 
888  chose = G4RandFlat::shoot(0.,2.*zArea+cArea);
889 
890  if( (chose>=0) && (chose < cArea) )
891  {
892  return G4ThreeVector (xRand,yRand,zRand);
893  }
894  else if( (chose >= cArea) && (chose < cArea + zArea) )
895  {
896  xRand = G4RandFlat::shoot(-1.*dx,dx);
897  yRand = std::sqrt(1.-sqr(xRand/dx));
898  yRand = G4RandFlat::shoot(-1.*yRand, yRand);
899  return G4ThreeVector (xRand,yRand,dz);
900  }
901  else
902  {
903  xRand = G4RandFlat::shoot(-1.*dx,dx);
904  yRand = std::sqrt(1.-sqr(xRand/dx));
905  yRand = G4RandFlat::shoot(-1.*yRand, yRand);
906  return G4ThreeVector (xRand,yRand,-1.*dz);
907  }
908 }
909 
910 
911 //
912 // CreatePolyhedron
913 //
915 {
916  // create cylinder with radius=1...
917  //
918  G4Polyhedron* eTube = new G4PolyhedronTube(0.,1.,dz);
919 
920  // apply non-uniform scaling...
921  //
922  eTube->Transform(G4Scale3D(dx,dy,1.));
923  return eTube;
924 }
925 
926 
927 //
928 // GetPolyhedron
929 //
931 {
932  if (!fpPolyhedron ||
933  fRebuildPolyhedron ||
935  fpPolyhedron->GetNumberOfRotationSteps())
936  {
937  G4AutoLock l(&polyhedronMutex);
938  delete fpPolyhedron;
939  fpPolyhedron = CreatePolyhedron();
940  fRebuildPolyhedron = false;
941  l.unlock();
942  }
943  return fpPolyhedron;
944 }
945 
946 
947 //
948 // DescribeYourselfTo
949 //
951 {
952  scene.AddSolid (*this);
953 }
954 
955 
956 //
957 // GetExtent
958 //
960 {
961  return G4VisExtent( -dx, dx, -dy, dy, -dz, dz );
962 }
void set(double x, double y, double z)
G4String GetName() const
const XML_Char * name
Definition: expat.h:151
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
ThreeVector shoot(const G4int Ap, const G4int Af)
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double x() const
G4double GetSurfaceArea()
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
virtual ~G4EllipticalTube()
const char * p
Definition: xmltok.h:285
virtual G4double GetCubicVolume()
Definition: G4VSolid.cc:189
G4ThreeVector GetPointOnSurface() const
G4VisExtent GetExtent() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
HepPolyhedron & Transform(const G4Transform3D &t)
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
double z() const
void DumpInfo() const
G4Polyhedron * GetPolyhedron() const
static constexpr double twopi
Definition: G4SIunits.hh:76
G4EllipticalTube & operator=(const G4EllipticalTube &rhs)
G4double CheckXY(const G4double x, const G4double y, const G4double toler) const
tuple b
Definition: test.py:12
G4int IntersectXY(const G4ThreeVector &p, const G4ThreeVector &v, G4double s[2]) const
G4double GetCubicVolume()
bool G4bool
Definition: G4Types.hh:79
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
HepGeom::Scale3D G4Scale3D
const G4int n
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
std::ostream & StreamInfo(std::ostream &os) const
std::vector< G4ThreeVector > G4ThreeVectorList
tuple v
Definition: test.py:18
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4Polyhedron * CreatePolyhedron() const
G4int G4Mutex
Definition: G4Threading.hh:173
static G4int GetNumberOfRotationSteps()
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Hep3Vector unit() const
double y() const
#define DBL_MIN
Definition: templates.hh:75
#define G4endl
Definition: G4ios.hh:61
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:111
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
static constexpr double pi
Definition: G4SIunits.hh:75
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4EllipticalTube(const G4String &name, G4double theDx, G4double theDy, G4double theDz)
#define DBL_MAX
Definition: templates.hh:83
G4GeometryType GetEntityType() const
virtual G4double GetSurfaceArea()
Definition: G4VSolid.cc:251
G4VSolid * Clone() const
EInside Inside(const G4ThreeVector &p) const