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