Geant4  10.02.p01
G4Ellipsoid.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 // $Id: G4Ellipsoid.cc 83572 2014-09-01 15:23:27Z gcosmo $
27 //
28 // class G4Ellipsoid
29 //
30 // Implementation for G4Ellipsoid class
31 //
32 // History:
33 //
34 // 10.11.99 G.Horton-Smith -- first writing, based on G4Sphere class
35 // 25.02.05 G.Guerrieri -- Modified for future Geant4 release
36 //
37 // --------------------------------------------------------------------
38 
39 #include "globals.hh"
40 
41 #include "G4Ellipsoid.hh"
42 
43 #include "G4VoxelLimits.hh"
44 #include "G4AffineTransform.hh"
45 #include "G4GeometryTolerance.hh"
46 
47 #include "meshdefs.hh"
48 #include "Randomize.hh"
49 
50 #include "G4VPVParameterisation.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 
65 //
66 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
67 // - note if pDPhi>2PI then reset to 2PI
68 
70  G4double pxSemiAxis,
71  G4double pySemiAxis,
72  G4double pzSemiAxis,
73  G4double pzBottomCut,
74  G4double pzTopCut)
75  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
76  fCubicVolume(0.), fSurfaceArea(0.), zBottomCut(0.), zTopCut(0.)
77 {
78  // note: for users that want to use the full ellipsoid it is useful
79  // to include a default for the cuts
80 
82 
85 
86  // Check Semi-Axis
87  if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
88  {
89  std::ostringstream message;
90  message << "Invalid semi-axis - " << GetName();
91  G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
92  FatalErrorInArgument, message);
93  }
94  SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis);
95 
96  if ( pzBottomCut == 0 && pzTopCut == 0 )
97  {
98  SetZCuts(-pzSemiAxis, pzSemiAxis);
99  }
100  else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
101  && (pzBottomCut < pzTopCut) )
102  {
103  SetZCuts(pzBottomCut, pzTopCut);
104  }
105  else
106  {
107  std::ostringstream message;
108  message << "Invalid z-coordinate for cutting plane - " << GetName();
109  G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
110  FatalErrorInArgument, message);
111  }
112 }
113 
115 //
116 // Fake default constructor - sets only member data and allocates memory
117 // for usage restricted to object persistency.
118 //
120  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0), kRadTolerance(0.),
121  halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
122  fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
123  semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
124 {
125 }
126 
128 //
129 // Destructor
130 
132 {
133  delete fpPolyhedron; fpPolyhedron = 0;
134 }
135 
137 //
138 // Copy constructor
139 
141  : G4VSolid(rhs),
142  fRebuildPolyhedron(false), fpPolyhedron(0),
143  kRadTolerance(rhs.kRadTolerance),
144  halfCarTolerance(rhs.halfCarTolerance),
145  halfRadTolerance(rhs.halfRadTolerance),
146  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
147  xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
148  zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
149  zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
150 {
151 }
152 
154 //
155 // Assignment operator
156 
158 {
159  // Check assignment to self
160  //
161  if (this == &rhs) { return *this; }
162 
163  // Copy base class data
164  //
165  G4VSolid::operator=(rhs);
166 
167  // Copy data
168  //
173  xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
175  zBottomCut = rhs.zBottomCut; zTopCut = rhs.zTopCut;
176  fRebuildPolyhedron = false;
177  delete fpPolyhedron; fpPolyhedron = 0;
178 
179  return *this;
180 }
181 
183 //
184 // Dispatch to parameterisation for replication mechanism dimension
185 // computation & modification.
186 
188  const G4int n,
189  const G4VPhysicalVolume* pRep)
190 {
191  p->ComputeDimensions(*this,n,pRep);
192 }
193 
195 //
196 // Calculate extent under transform and specified limit
197 
198 G4bool
200  const G4VoxelLimits& pVoxelLimit,
201  const G4AffineTransform& pTransform,
202  G4double& pMin, G4double& pMax) const
203 {
204  if (!pTransform.IsRotated())
205  {
206  // Special case handling for unrotated solid ellipsoid
207  // Compute x/y/z mins and maxs for bounding box respecting limits,
208  // with early returns if outside limits. Then switch() on pAxis,
209  // and compute exact x and y limit for x/y case
210 
211  G4double xoffset,xMin,xMax;
212  G4double yoffset,yMin,yMax;
213  G4double zoffset,zMin,zMax;
214 
215  G4double maxDiff,newMin,newMax;
216  G4double xoff,yoff;
217 
218  xoffset=pTransform.NetTranslation().x();
219  xMin=xoffset - xSemiAxis;
220  xMax=xoffset + xSemiAxis;
221  if (pVoxelLimit.IsXLimited())
222  {
223  if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
224  || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
225  {
226  return false;
227  }
228  else
229  {
230  if (xMin<pVoxelLimit.GetMinXExtent())
231  {
232  xMin=pVoxelLimit.GetMinXExtent();
233  }
234  if (xMax>pVoxelLimit.GetMaxXExtent())
235  {
236  xMax=pVoxelLimit.GetMaxXExtent();
237  }
238  }
239  }
240 
241  yoffset=pTransform.NetTranslation().y();
242  yMin=yoffset - ySemiAxis;
243  yMax=yoffset + ySemiAxis;
244  if (pVoxelLimit.IsYLimited())
245  {
246  if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
247  || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
248  {
249  return false;
250  }
251  else
252  {
253  if (yMin<pVoxelLimit.GetMinYExtent())
254  {
255  yMin=pVoxelLimit.GetMinYExtent();
256  }
257  if (yMax>pVoxelLimit.GetMaxYExtent())
258  {
259  yMax=pVoxelLimit.GetMaxYExtent();
260  }
261  }
262  }
263 
264  zoffset=pTransform.NetTranslation().z();
265  zMin=zoffset + (-zSemiAxis > zBottomCut ? -zSemiAxis : zBottomCut);
266  zMax=zoffset + ( zSemiAxis < zTopCut ? zSemiAxis : zTopCut);
267  if (pVoxelLimit.IsZLimited())
268  {
269  if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
270  || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
271  {
272  return false;
273  }
274  else
275  {
276  if (zMin<pVoxelLimit.GetMinZExtent())
277  {
278  zMin=pVoxelLimit.GetMinZExtent();
279  }
280  if (zMax>pVoxelLimit.GetMaxZExtent())
281  {
282  zMax=pVoxelLimit.GetMaxZExtent();
283  }
284  }
285  }
286 
287  // if here, then known to cut bounding box around ellipsoid
288  //
289  xoff = (xoffset < xMin) ? (xMin-xoffset)
290  : (xoffset > xMax) ? (xoffset-xMax) : 0.0;
291  yoff = (yoffset < yMin) ? (yMin-yoffset)
292  : (yoffset > yMax) ? (yoffset-yMax) : 0.0;
293 
294  // detailed calculations
295  // NOTE: does not use X or Y offsets to adjust Z range,
296  // and does not use Z offset to adjust X or Y range,
297  // which is consistent with G4Sphere::CalculateExtent behavior
298  //
299  switch (pAxis)
300  {
301  case kXAxis:
302  if (yoff==0.)
303  {
304  // YZ limits cross max/min x => no change
305  //
306  pMin=xMin;
307  pMax=xMax;
308  }
309  else
310  {
311  // YZ limits don't cross max/min x => compute max delta x,
312  // hence new mins/maxs
313  //
314  maxDiff= 1.0-sqr(yoff/ySemiAxis);
315  if (maxDiff < 0.0) { return false; }
316  maxDiff= xSemiAxis * std::sqrt(maxDiff);
317  newMin=xoffset-maxDiff;
318  newMax=xoffset+maxDiff;
319  pMin=(newMin<xMin) ? xMin : newMin;
320  pMax=(newMax>xMax) ? xMax : newMax;
321  }
322  break;
323  case kYAxis:
324  if (xoff==0.)
325  {
326  // XZ limits cross max/min y => no change
327  //
328  pMin=yMin;
329  pMax=yMax;
330  }
331  else
332  {
333  // XZ limits don't cross max/min y => compute max delta y,
334  // hence new mins/maxs
335  //
336  maxDiff= 1.0-sqr(xoff/xSemiAxis);
337  if (maxDiff < 0.0) { return false; }
338  maxDiff= ySemiAxis * std::sqrt(maxDiff);
339  newMin=yoffset-maxDiff;
340  newMax=yoffset+maxDiff;
341  pMin=(newMin<yMin) ? yMin : newMin;
342  pMax=(newMax>yMax) ? yMax : newMax;
343  }
344  break;
345  case kZAxis:
346  pMin=zMin;
347  pMax=zMax;
348  break;
349  default:
350  break;
351  }
352 
353  pMin-=kCarTolerance;
354  pMax+=kCarTolerance;
355  return true;
356  }
357  else // not rotated
358  {
359  G4int i,j,noEntries,noBetweenSections;
360  G4bool existsAfterClip=false;
361 
362  // Calculate rotated vertex coordinates
363 
364  G4int noPolygonVertices=0;
365  G4ThreeVectorList* vertices =
366  CreateRotatedVertices(pTransform,noPolygonVertices);
367 
368  pMin=+kInfinity;
369  pMax=-kInfinity;
370 
371  noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections
372  noBetweenSections=noEntries-noPolygonVertices;
373 
374  G4ThreeVectorList ThetaPolygon;
375  for (i=0;i<noEntries;i+=noPolygonVertices)
376  {
377  for(j=0;j<(noPolygonVertices/2)-1;j++)
378  {
379  ThetaPolygon.push_back((*vertices)[i+j]);
380  ThetaPolygon.push_back((*vertices)[i+j+1]);
381  ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]);
382  ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]);
383  CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
384  ThetaPolygon.clear();
385  }
386  }
387  for (i=0;i<noBetweenSections;i+=noPolygonVertices)
388  {
389  for(j=0;j<noPolygonVertices-1;j++)
390  {
391  ThetaPolygon.push_back((*vertices)[i+j]);
392  ThetaPolygon.push_back((*vertices)[i+j+1]);
393  ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]);
394  ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]);
395  CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
396  ThetaPolygon.clear();
397  }
398  ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]);
399  ThetaPolygon.push_back((*vertices)[i]);
400  ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]);
401  ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]);
402  CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
403  ThetaPolygon.clear();
404  }
405  if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
406  {
407  existsAfterClip=true;
408 
409  // Add 2*tolerance to avoid precision troubles
410  //
411  pMin-=kCarTolerance;
412  pMax+=kCarTolerance;
413 
414  }
415  else
416  {
417  // Check for case where completely enveloping clipping volume
418  // If point inside then we are confident that the solid completely
419  // envelopes the clipping volume. Hence set min/max extents according
420  // to clipping volume extents along the specified axis.
421  //
423  clipCentre((pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
424  (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
425  (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
426 
427  if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
428  {
429  existsAfterClip=true;
430  pMin=pVoxelLimit.GetMinExtent(pAxis);
431  pMax=pVoxelLimit.GetMaxExtent(pAxis);
432  }
433  }
434  delete vertices;
435  return existsAfterClip;
436  }
437 }
438 
440 //
441 // Return whether point inside/outside/on surface
442 // Split into radius, phi, theta checks
443 // Each check modifies `in', or returns as approprate
444 
446 {
447  G4double rad2oo, // outside surface outer tolerance
448  rad2oi; // outside surface inner tolerance
449  EInside in;
450 
451  // check this side of z cut first, because that's fast
452  //
453  if (p.z() < zBottomCut-halfRadTolerance) { return in=kOutside; }
454  if (p.z() > zTopCut+halfRadTolerance) { return in=kOutside; }
455 
456  rad2oo= sqr(p.x()/(xSemiAxis+halfRadTolerance))
457  + sqr(p.y()/(ySemiAxis+halfRadTolerance))
458  + sqr(p.z()/(zSemiAxis+halfRadTolerance));
459 
460  if (rad2oo > 1.0) { return in=kOutside; }
461 
462  rad2oi= sqr(p.x()*(1.0+halfRadTolerance/xSemiAxis)/xSemiAxis)
463  + sqr(p.y()*(1.0+halfRadTolerance/ySemiAxis)/ySemiAxis)
464  + sqr(p.z()*(1.0+halfRadTolerance/zSemiAxis)/zSemiAxis);
465 
466  // Check radial surfaces
467  // sets `in' (already checked for rad2oo > 1.0)
468  //
469  if (rad2oi < 1.0)
470  {
471  in = ( (p.z() < zBottomCut+halfRadTolerance)
472  || (p.z() > zTopCut-halfRadTolerance) ) ? kSurface : kInside;
473  if ( rad2oi > 1.0-halfRadTolerance ) { in=kSurface; }
474  }
475  else
476  {
477  in = kSurface;
478  }
479  return in;
480 
481 }
482 
484 //
485 // Return unit normal of surface closest to p not protected against p=0
486 
488 {
489  G4double distR, distZBottom, distZTop;
490 
491  // normal vector with special magnitude: parallel to normal, units 1/length
492  // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside
493  //
494  G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
495  p.y()/(ySemiAxis*ySemiAxis),
496  p.z()/(zSemiAxis*zSemiAxis));
497  G4double radius = 1.0/norm.mag();
498 
499  // approximate distance to curved surface
500  //
501  distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
502 
503  // Distance to z-cut plane
504  //
505  distZBottom = std::fabs( p.z() - zBottomCut );
506  distZTop = std::fabs( p.z() - zTopCut );
507 
508  if ( (distZBottom < distR) || (distZTop < distR) )
509  {
510  return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
511  }
512  return ( norm *= radius );
513 }
514 
516 //
517 // Calculate distance to shape from outside, along normalised vector
518 // - return kInfinity if no intersection, or intersection distance <= tolerance
519 //
520 
522  const G4ThreeVector& v ) const
523 {
525  const G4double dRmax = 100.*std::min(distMin,zSemiAxis);
526  distMin= kInfinity;
527 
528  // check to see if Z plane is relevant
529  if (p.z() <= zBottomCut+halfCarTolerance)
530  {
531  if (v.z() <= 0.0) { return distMin; }
532  G4double distZ = (zBottomCut - p.z()) / v.z();
533 
534  if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
535  {
536  // early exit since can't intercept curved surface if we reach here
537  if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
538  return distMin= distZ;
539  }
540  }
541  if (p.z() >= zTopCut-halfCarTolerance)
542  {
543  if (v.z() >= 0.0) { return distMin;}
544  G4double distZ = (zTopCut - p.z()) / v.z();
545  if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
546  {
547  // early exit since can't intercept curved surface if we reach here
548  if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
549  return distMin= distZ;
550  }
551  }
552  // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface
553 
554  // now check curved surface intercept
555  G4double A,B,C;
556 
557  A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
558  C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0;
559  B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis)
560  + p.y()*v.y()/(ySemiAxis*ySemiAxis)
561  + p.z()*v.z()/(zSemiAxis*zSemiAxis) );
562 
563  C= B*B - 4.0*A*C;
564  if (C > 0.0)
565  {
566  G4double distR= (-B - std::sqrt(C)) / (2.0*A);
567  G4double intZ = p.z()+distR*v.z();
568  if ( (distR > halfRadTolerance)
569  && (intZ >= zBottomCut-halfRadTolerance)
570  && (intZ <= zTopCut+halfRadTolerance) )
571  {
572  distMin = distR;
573  }
574  else if( (distR >- halfRadTolerance)
575  && (intZ >= zBottomCut-halfRadTolerance)
576  && (intZ <= zTopCut+halfRadTolerance) )
577  {
578  // p is on the curved surface, DistanceToIn returns 0 or kInfinity:
579  // DistanceToIn returns 0, if second root is positive (means going inside)
580  // If second root is negative, DistanceToIn returns kInfinity (outside)
581  //
582  distR = (-B + std::sqrt(C) ) / (2.0*A);
583  if(distR>0.) { distMin=0.; }
584  }
585  else
586  {
587  distR= (-B + std::sqrt(C)) / (2.0*A);
588  intZ = p.z()+distR*v.z();
589  if ( (distR > halfRadTolerance)
590  && (intZ >= zBottomCut-halfRadTolerance)
591  && (intZ <= zTopCut+halfRadTolerance) )
592  {
594  if (norm.dot(v)<0.) { distMin = distR; }
595  }
596  }
597  if ( (distMin!=kInfinity) && (distMin>dRmax) )
598  { // Avoid rounding errors due to precision issues on
599  // 64 bits systems. Split long distances and recompute
600  G4double fTerm = distMin-std::fmod(distMin,dRmax);
601  distMin = fTerm + DistanceToIn(p+fTerm*v,v);
602  }
603  }
604 
605  if (std::fabs(distMin)<halfRadTolerance) { distMin=0.; }
606  return distMin;
607 }
608 
610 //
611 // Calculate distance (<= actual) to closest surface of shape from outside
612 // - Return 0 if point inside
613 
615 {
616  G4double distR, distZ;
617 
618  // normal vector: parallel to normal, magnitude 1/(characteristic radius)
619  //
620  G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
621  p.y()/(ySemiAxis*ySemiAxis),
622  p.z()/(zSemiAxis*zSemiAxis));
623  G4double radius= 1.0/norm.mag();
624 
625  // approximate distance to curved surface ( <= actual distance )
626  //
627  distR= (p*norm - 1.0) * radius / 2.0;
628 
629  // Distance to z-cut plane
630  //
631  distZ= zBottomCut - p.z();
632  if (distZ < 0.0)
633  {
634  distZ = p.z() - zTopCut;
635  }
636 
637  // Distance to closest surface from outside
638  //
639  if (distZ < 0.0)
640  {
641  return (distR < 0.0) ? 0.0 : distR;
642  }
643  else if (distR < 0.0)
644  {
645  return distZ;
646  }
647  else
648  {
649  return (distZ < distR) ? distZ : distR;
650  }
651 }
652 
654 //
655 // Calculate distance to surface of shape from `inside', allowing for tolerance
656 
658  const G4ThreeVector& v,
659  const G4bool calcNorm,
660  G4bool *validNorm,
661  G4ThreeVector *n ) const
662 {
663  G4double distMin;
664  enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
665 
666  distMin= kInfinity;
667  surface= kNoSurf;
668 
669  // check to see if Z plane is relevant
670  //
671  if (v.z() < 0.0)
672  {
673  G4double distZ = (zBottomCut - p.z()) / v.z();
674  if (distZ < 0.0)
675  {
676  distZ= 0.0;
677  if (!calcNorm) {return 0.0;}
678  }
679  distMin= distZ;
680  surface= kPlaneSurf;
681  }
682  if (v.z() > 0.0)
683  {
684  G4double distZ = (zTopCut - p.z()) / v.z();
685  if (distZ < 0.0)
686  {
687  distZ= 0.0;
688  if (!calcNorm) {return 0.0;}
689  }
690  distMin= distZ;
691  surface= kPlaneSurf;
692  }
693 
694  // normal vector: parallel to normal, magnitude 1/(characteristic radius)
695  //
696  G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis),
697  p.y()/(ySemiAxis*ySemiAxis),
698  p.z()/(zSemiAxis*zSemiAxis));
699 
700  // now check curved surface intercept
701  //
702  G4double A,B,C;
703 
704  A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
705  C= (p * nearnorm) - 1.0;
706  B= 2.0 * (v * nearnorm);
707 
708  C= B*B - 4.0*A*C;
709  if (C > 0.0)
710  {
711  G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
712  if (distR < 0.0)
713  {
714  distR= 0.0;
715  if (!calcNorm) {return 0.0;}
716  }
717  if (distR < distMin)
718  {
719  distMin= distR;
720  surface= kCurvedSurf;
721  }
722  }
723 
724  // set normal if requested
725  //
726  if (calcNorm)
727  {
728  if (surface == kNoSurf)
729  {
730  *validNorm = false;
731  }
732  else
733  {
734  *validNorm = true;
735  switch (surface)
736  {
737  case kPlaneSurf:
738  *n= G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
739  break;
740  case kCurvedSurf:
741  {
742  G4ThreeVector pexit= p + distMin*v;
743  G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis),
744  pexit.y()/(ySemiAxis*ySemiAxis),
745  pexit.z()/(zSemiAxis*zSemiAxis));
746  truenorm *= 1.0/truenorm.mag();
747  *n= truenorm;
748  } break;
749  default: // Should never reach this case ...
750  DumpInfo();
751  std::ostringstream message;
752  G4int oldprc = message.precision(16);
753  message << "Undefined side for valid surface normal to solid."
754  << G4endl
755  << "Position:" << G4endl
756  << " p.x() = " << p.x()/mm << " mm" << G4endl
757  << " p.y() = " << p.y()/mm << " mm" << G4endl
758  << " p.z() = " << p.z()/mm << " mm" << G4endl
759  << "Direction:" << G4endl << G4endl
760  << " v.x() = " << v.x() << G4endl
761  << " v.y() = " << v.y() << G4endl
762  << " v.z() = " << v.z() << G4endl
763  << "Proposed distance :" << G4endl
764  << " distMin = " << distMin/mm << " mm";
765  message.precision(oldprc);
766  G4Exception("G4Ellipsoid::DistanceToOut(p,v,..)",
767  "GeomSolids1002", JustWarning, message);
768  break;
769  }
770  }
771  }
772 
773  return distMin;
774 }
775 
777 //
778 // Calculate distance (<=actual) to closest surface of shape from inside
779 
781 {
782  G4double distR, distZ;
783 
784 #ifdef G4SPECSDEBUG
785  if( Inside(p) == kOutside )
786  {
787  DumpInfo();
788  std::ostringstream message;
789  G4int oldprc = message.precision(16);
790  message << "Point p is outside !?" << G4endl
791  << "Position:" << G4endl
792  << " p.x() = " << p.x()/mm << " mm" << G4endl
793  << " p.y() = " << p.y()/mm << " mm" << G4endl
794  << " p.z() = " << p.z()/mm << " mm";
795  message.precision(oldprc) ;
796  G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
797  JustWarning, message);
798  }
799 #endif
800 
801  // Normal vector: parallel to normal, magnitude 1/(characteristic radius)
802  //
803  G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
804  p.y()/(ySemiAxis*ySemiAxis),
805  p.z()/(zSemiAxis*zSemiAxis));
806 
807  // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag())
808  //
809  G4double radius= p.mag();
810  G4double tmp= norm.mag();
811  if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
812 
813  // Approximate distance to curved surface ( <= actual distance )
814  //
815  distR = (1.0 - p*norm) * radius / 2.0;
816 
817  // Distance to z-cut plane
818  //
819  distZ = p.z() - zBottomCut;
820  if (distZ < 0.0) {distZ= zTopCut - p.z();}
821 
822  // Distance to closest surface from inside
823  //
824  if ( (distZ < 0.0) || (distR < 0.0) )
825  {
826  return 0.0;
827  }
828  else
829  {
830  return (distZ < distR) ? distZ : distR;
831  }
832 }
833 
835 //
836 // Create a List containing the transformed vertices
837 // Ordering [0-3] -fDz cross section
838 // [4-7] +fDz cross section such that [0] is below [4],
839 // [1] below [5] etc.
840 // Note:
841 // Caller has deletion resposibility
842 // Potential improvement: For last slice, use actual ending angle
843 // to avoid rounding error problems.
844 
847  G4int& noPolygonVertices) const
848 {
849  G4ThreeVectorList *vertices;
850  G4ThreeVector vertex;
851  G4double meshAnglePhi, meshRMaxFactor,
852  crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
853  G4double meshTheta, crossTheta, startTheta;
854  G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
855  G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
856 
857  // Phi cross sections
858  //
859  noPhiCrossSections=G4int (twopi/kMeshAngleDefault)+1; // = 9!
860 
861 /*
862  if (noPhiCrossSections<kMinMeshSections) // <3
863  {
864  noPhiCrossSections=kMinMeshSections;
865  }
866  else if (noPhiCrossSections>kMaxMeshSections) // >37
867  {
868  noPhiCrossSections=kMaxMeshSections;
869  }
870 */
871  meshAnglePhi=twopi/(noPhiCrossSections-1);
872 
873  // Set start angle such that mesh will be at fRMax
874  // on the x axis. Will give better extent calculations when not rotated.
875 
876  sAnglePhi = -meshAnglePhi*0.5;
877 
878  // Theta cross sections
879 
880  noThetaSections = G4int(pi/kMeshAngleDefault)+3; // = 7!
881 
882 /*
883  if (noThetaSections<kMinMeshSections) // <3
884  {
885  noThetaSections=kMinMeshSections;
886  }
887  else if (noThetaSections>kMaxMeshSections) // >37
888  {
889  noThetaSections=kMaxMeshSections;
890  }
891 */
892  meshTheta= pi/(noThetaSections-2);
893 
894  // Set start angle such that mesh will be at fRMax
895  // on the z axis. Will give better extent calculations when not rotated.
896 
897  startTheta = -meshTheta*0.5;
898 
899  meshRMaxFactor = 1.0/std::cos(0.5*
900  std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
901  rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis);
902  if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis;
903  rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
904  rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
905  rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
906  G4double* cosCrossTheta = new G4double[noThetaSections];
907  G4double* sinCrossTheta = new G4double[noThetaSections];
908  vertices=new G4ThreeVectorList(noPhiCrossSections*noThetaSections);
909  if (vertices && cosCrossTheta && sinCrossTheta)
910  {
911  for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
912  crossSectionTheta++)
913  {
914  // Compute sine and cosine table (for historical reasons)
915  //
916  crossTheta=startTheta+crossSectionTheta*meshTheta;
917  cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
918  sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
919  }
920  for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
921  crossSectionPhi++)
922  {
923  crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
924  coscrossAnglePhi=std::cos(crossAnglePhi);
925  sincrossAnglePhi=std::sin(crossAnglePhi);
926  for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
927  crossSectionTheta++)
928  {
929  // Compute coordinates of cross section at section crossSectionPhi
930  //
931  rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
932  ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
933  rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
934  if (rz < zBottomCut)
935  { rz= zBottomCut; }
936  if (rz > zTopCut)
937  { rz= zTopCut; }
938  vertex= G4ThreeVector(rx,ry,rz);
939  vertices->push_back(pTransform.TransformPoint(vertex));
940  } // Theta forward
941  } // Phi
942  noPolygonVertices = noThetaSections ;
943  }
944  else
945  {
946  DumpInfo();
947  G4Exception("G4Ellipsoid::CreateRotatedVertices()",
948  "GeomSolids0003", FatalException,
949  "Error in allocation of vertices. Out of memory !");
950  }
951 
952  delete[] cosCrossTheta;
953  delete[] sinCrossTheta;
954 
955  return vertices;
956 }
957 
959 //
960 // G4EntityType
961 
963 {
964  return G4String("G4Ellipsoid");
965 }
966 
968 //
969 // Make a clone of the object
970 
972 {
973  return new G4Ellipsoid(*this);
974 }
975 
977 //
978 // Stream object contents to an output stream
979 
980 std::ostream& G4Ellipsoid::StreamInfo( std::ostream& os ) const
981 {
982  G4int oldprc = os.precision(16);
983  os << "-----------------------------------------------------------\n"
984  << " *** Dump for solid - " << GetName() << " ***\n"
985  << " ===================================================\n"
986  << " Solid type: G4Ellipsoid\n"
987  << " Parameters: \n"
988 
989  << " semi-axis x: " << xSemiAxis/mm << " mm \n"
990  << " semi-axis y: " << ySemiAxis/mm << " mm \n"
991  << " semi-axis z: " << zSemiAxis/mm << " mm \n"
992  << " max semi-axis: " << semiAxisMax/mm << " mm \n"
993  << " lower cut plane level z: " << zBottomCut/mm << " mm \n"
994  << " upper cut plane level z: " << zTopCut/mm << " mm \n"
995  << "-----------------------------------------------------------\n";
996  os.precision(oldprc);
997 
998  return os;
999 }
1000 
1002 //
1003 // GetPointOnSurface
1004 
1006 {
1007  G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
1008  G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3;
1009 
1010  max1 = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
1011  max1 = max1 > zSemiAxis ? max1 : zSemiAxis;
1012  if (max1 == xSemiAxis) { max2 = ySemiAxis; max3 = zSemiAxis; }
1013  else if (max1 == ySemiAxis) { max2 = xSemiAxis; max3 = zSemiAxis; }
1014  else { max2 = xSemiAxis; max3 = ySemiAxis; }
1015 
1016  phi = RandFlat::shoot(0.,twopi);
1017 
1018  cosphi = std::cos(phi); sinphi = std::sin(phi);
1020  sintheta = std::sqrt(1.-sqr(costheta));
1021 
1022  alpha = 1.-sqr(max2/max1); beta = 1.-sqr(max3/max1);
1023 
1024  aTop = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis));
1025  aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis));
1026 
1027  // approximation
1028  // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf"
1029  aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)-
1030  1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta)));
1031 
1032  aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis);
1033 
1034  if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis )
1035  || ( zTopCut == 0 && zBottomCut ==0 ) )
1036  {
1037  aTop = 0; aBottom = 0;
1038  }
1039 
1040  chose = RandFlat::shoot(0.,aTop + aBottom + aCurved);
1041 
1042  if(chose < aCurved)
1043  {
1044  xRand = xSemiAxis*sintheta*cosphi;
1045  yRand = ySemiAxis*sintheta*sinphi;
1046  zRand = zSemiAxis*costheta;
1047  return G4ThreeVector (xRand,yRand,zRand);
1048  }
1049  else if(chose >= aCurved && chose < aCurved + aTop)
1050  {
1051  xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
1052  * std::sqrt(1-sqr(zTopCut/zSemiAxis));
1053  yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
1054  * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis));
1055  zRand = zTopCut;
1056  return G4ThreeVector (xRand,yRand,zRand);
1057  }
1058  else
1059  {
1060  xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
1061  * std::sqrt(1-sqr(zBottomCut/zSemiAxis));
1062  yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
1063  * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis));
1064  zRand = zBottomCut;
1065  return G4ThreeVector (xRand,yRand,zRand);
1066  }
1067 }
1068 
1070 //
1071 // Methods for visualisation
1072 
1074 {
1075  scene.AddSolid(*this);
1076 }
1077 
1079 {
1080  // Define the sides of the box into which the G4Ellipsoid instance would fit.
1081  //
1085 }
1086 
1088 {
1090  zBottomCut, zTopCut);
1091 }
1092 
1094 {
1095  if (!fpPolyhedron ||
1098  fpPolyhedron->GetNumberOfRotationSteps())
1099  {
1100  G4AutoLock l(&polyhedronMutex);
1101  delete fpPolyhedron;
1103  fRebuildPolyhedron = false;
1104  l.unlock();
1105  }
1106  return fpPolyhedron;
1107 }
G4Polyhedron * GetPolyhedron() const
G4double zSemiAxis
Definition: G4Ellipsoid.hh:146
G4String GetName() const
G4double zTopCut
Definition: G4Ellipsoid.hh:146
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double halfRadTolerance
Definition: G4Ellipsoid.hh:142
G4double fSurfaceArea
Definition: G4Ellipsoid.hh:145
G4GeometryType GetEntityType() const
Definition: G4Ellipsoid.cc:962
G4bool fRebuildPolyhedron
Definition: G4Ellipsoid.hh:136
G4double xSemiAxis
Definition: G4Ellipsoid.hh:146
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4AffineTransform Inverse() const
G4bool IsYLimited() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Ellipsoid.cc:657
EInside Inside(const G4ThreeVector &p) const
Definition: G4Ellipsoid.cc:445
G4bool IsRotated() const
G4ThreeVector NetTranslation() const
G4bool IsXLimited() const
G4double a
Definition: TRTMaterials.hh:39
double C(double temp)
double B(double temperature)
virtual ~G4Ellipsoid()
Definition: G4Ellipsoid.cc:131
virtual void AddSolid(const G4Box &)=0
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Ellipsoid.cc:199
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
void DumpInfo() const
G4double ySemiAxis
Definition: G4Ellipsoid.hh:146
G4double kRadTolerance
Definition: G4Ellipsoid.hh:141
void SetZCuts(G4double newzBottomCut, G4double newzTopCut)
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4Ellipsoid(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double pzSemiAxis, G4double pzBottomCut=0, G4double pzTopCut=0)
Definition: G4Ellipsoid.cc:69
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
double A(double temperature)
void CalculateClippedPolygonExtent(G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:427
G4double semiAxisMax
Definition: G4Ellipsoid.hh:146
G4VisExtent GetExtent() const
G4Polyhedron * fpPolyhedron
Definition: G4Ellipsoid.hh:137
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pT, G4int &noPV) const
Definition: G4Ellipsoid.cc:846
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
G4Ellipsoid & operator=(const G4Ellipsoid &rhs)
Definition: G4Ellipsoid.cc:157
G4Polyhedron * CreatePolyhedron() const
G4VSolid * Clone() const
Definition: G4Ellipsoid.cc:971
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Ellipsoid.cc:521
G4double zBottomCut
Definition: G4Ellipsoid.hh:146
G4double fCubicVolume
Definition: G4Ellipsoid.hh:144
const G4int n
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Ellipsoid.cc:487
G4double halfCarTolerance
Definition: G4Ellipsoid.hh:142
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double GetMinXExtent() const
G4int G4Mutex
Definition: G4Threading.hh:173
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Ellipsoid.cc:980
G4double GetMaxZExtent() const
static const double pi
Definition: G4SIunits.hh:74
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetSemiAxis(G4double x, G4double y, G4double z)
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
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
G4ThreeVector GetPointOnSurface() const
double G4double
Definition: G4Types.hh:76
G4double GetMaxExtent(const EAxis pAxis) const
static const G4double alpha
G4bool IsZLimited() const
static const double mm
Definition: G4SIunits.hh:114
G4double GetMinExtent(const EAxis pAxis) const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Ellipsoid.cc:187
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42
static G4GeometryTolerance * GetInstance()
static int max3(int a, int b, int c)