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