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