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