Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 100819 2016-11-02 15:17:36Z 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 // 26.10.16 E.Tcherniaev: Added Extent(pmin,pmax),
37 // reimplemented CalculateExtent() using G4BoundingEnvelope,
38 // removed CreateRotatedVertices()
39 //
40 // --------------------------------------------------------------------
41 
42 #include "globals.hh"
43 
44 #include "G4Ellipsoid.hh"
45 
46 #include "G4VoxelLimits.hh"
47 #include "G4AffineTransform.hh"
48 #include "G4GeometryTolerance.hh"
49 #include "G4BoundingEnvelope.hh"
50 
51 #include "meshdefs.hh"
52 #include "Randomize.hh"
53 
54 #include "G4VPVParameterisation.hh"
55 
56 #include "G4VGraphicsScene.hh"
57 #include "G4VisExtent.hh"
58 
59 #include "G4AutoLock.hh"
60 
61 namespace
62 {
63  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
64 }
65 
66 using namespace CLHEP;
67 
69 //
70 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
71 // - note if pDPhi>2PI then reset to 2PI
72 
74  G4double pxSemiAxis,
75  G4double pySemiAxis,
76  G4double pzSemiAxis,
77  G4double pzBottomCut,
78  G4double pzTopCut)
79  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
80  fCubicVolume(0.), fSurfaceArea(0.), zBottomCut(0.), zTopCut(0.)
81 {
82  // note: for users that want to use the full ellipsoid it is useful
83  // to include a default for the cuts
84 
86 
87  halfCarTolerance = kCarTolerance*0.5;
88  halfRadTolerance = kRadTolerance*0.5;
89 
90  // Check Semi-Axis
91  if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
92  {
93  std::ostringstream message;
94  message << "Invalid semi-axis - " << GetName();
95  G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
96  FatalErrorInArgument, message);
97  }
98  SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis);
99 
100  if ( pzBottomCut == 0 && pzTopCut == 0 )
101  {
102  SetZCuts(-pzSemiAxis, pzSemiAxis);
103  }
104  else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
105  && (pzBottomCut < pzTopCut) )
106  {
107  SetZCuts(pzBottomCut, pzTopCut);
108  }
109  else
110  {
111  std::ostringstream message;
112  message << "Invalid z-coordinate for cutting plane - " << GetName();
113  G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
114  FatalErrorInArgument, message);
115  }
116 }
117 
119 //
120 // Fake default constructor - sets only member data and allocates memory
121 // for usage restricted to object persistency.
122 //
124  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0), kRadTolerance(0.),
125  halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
126  fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
127  semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
128 {
129 }
130 
132 //
133 // Destructor
134 
136 {
137  delete fpPolyhedron; fpPolyhedron = 0;
138 }
139 
141 //
142 // Copy constructor
143 
145  : G4VSolid(rhs),
146  fRebuildPolyhedron(false), fpPolyhedron(0),
147  kRadTolerance(rhs.kRadTolerance),
148  halfCarTolerance(rhs.halfCarTolerance),
149  halfRadTolerance(rhs.halfRadTolerance),
150  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
151  xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
152  zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
153  zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
154 {
155 }
156 
158 //
159 // Assignment operator
160 
162 {
163  // Check assignment to self
164  //
165  if (this == &rhs) { return *this; }
166 
167  // Copy base class data
168  //
169  G4VSolid::operator=(rhs);
170 
171  // Copy data
172  //
173  kRadTolerance = rhs.kRadTolerance;
174  halfCarTolerance = rhs.halfCarTolerance;
175  halfRadTolerance = rhs.halfRadTolerance;
176  fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
177  xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
178  zSemiAxis = rhs.zSemiAxis; semiAxisMax = rhs.semiAxisMax;
179  zBottomCut = rhs.zBottomCut; zTopCut = rhs.zTopCut;
180  fRebuildPolyhedron = false;
181  delete fpPolyhedron; fpPolyhedron = 0;
182 
183  return *this;
184 }
185 
187 //
188 // Dispatch to parameterisation for replication mechanism dimension
189 // computation & modification.
190 
192  const G4int n,
193  const G4VPhysicalVolume* pRep)
194 {
195  p->ComputeDimensions(*this,n,pRep);
196 }
197 
199 //
200 // Get bounding box
201 
203 {
204  G4double dx = GetSemiAxisMax(0);
205  G4double dy = GetSemiAxisMax(1);
206  G4double dz = GetSemiAxisMax(2);
207  G4double zmin = std::max(-dz,GetZBottomCut());
208  G4double zmax = std::min( dz,GetZTopCut());
209  pMin.set(-dx,-dy,zmin);
210  pMax.set( dx, dy,zmax);
211 
212  // Check correctness of the bounding box
213  //
214  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
215  {
216  std::ostringstream message;
217  message << "Bad bounding box (min >= max) for solid: "
218  << GetName() << " !"
219  << "\npMin = " << pMin
220  << "\npMax = " << pMax;
221  G4Exception("G4Ellipsoid::Extent()", "GeomMgt0001", JustWarning, message);
222  DumpInfo();
223  }
224 }
225 
227 //
228 // Calculate extent under transform and specified limit
229 
230 G4bool
232  const G4VoxelLimits& pVoxelLimit,
233  const G4AffineTransform& pTransform,
234  G4double& pMin, G4double& pMax) const
235 {
236  G4ThreeVector bmin, bmax;
237 
238  // Get bounding box
239  Extent(bmin,bmax);
240 
241  // Find extent
242  G4BoundingEnvelope bbox(bmin,bmax);
243  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
244 }
245 
247 //
248 // Return whether point inside/outside/on surface
249 // Split into radius, phi, theta checks
250 // Each check modifies `in', or returns as approprate
251 
253 {
254  G4double rad2oo, // outside surface outer tolerance
255  rad2oi; // outside surface inner tolerance
256  EInside in;
257 
258  // check this side of z cut first, because that's fast
259  //
260  if (p.z() < zBottomCut-halfRadTolerance) { return in=kOutside; }
261  if (p.z() > zTopCut+halfRadTolerance) { return in=kOutside; }
262 
263  rad2oo= sqr(p.x()/(xSemiAxis+halfRadTolerance))
264  + sqr(p.y()/(ySemiAxis+halfRadTolerance))
265  + sqr(p.z()/(zSemiAxis+halfRadTolerance));
266 
267  if (rad2oo > 1.0) { return in=kOutside; }
268 
269  rad2oi= sqr(p.x()*(1.0+halfRadTolerance/xSemiAxis)/xSemiAxis)
270  + sqr(p.y()*(1.0+halfRadTolerance/ySemiAxis)/ySemiAxis)
271  + sqr(p.z()*(1.0+halfRadTolerance/zSemiAxis)/zSemiAxis);
272 
273  // Check radial surfaces
274  // sets `in' (already checked for rad2oo > 1.0)
275  //
276  if (rad2oi < 1.0)
277  {
278  in = ( (p.z() < zBottomCut+halfRadTolerance)
279  || (p.z() > zTopCut-halfRadTolerance) ) ? kSurface : kInside;
280  if ( rad2oi > 1.0-halfRadTolerance ) { in=kSurface; }
281  }
282  else
283  {
284  in = kSurface;
285  }
286  return in;
287 
288 }
289 
291 //
292 // Return unit normal of surface closest to p not protected against p=0
293 
295 {
296  G4double distR, distZBottom, distZTop;
297 
298  // normal vector with special magnitude: parallel to normal, units 1/length
299  // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside
300  //
301  G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
302  p.y()/(ySemiAxis*ySemiAxis),
303  p.z()/(zSemiAxis*zSemiAxis));
304  G4double radius = 1.0/norm.mag();
305 
306  // approximate distance to curved surface
307  //
308  distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
309 
310  // Distance to z-cut plane
311  //
312  distZBottom = std::fabs( p.z() - zBottomCut );
313  distZTop = std::fabs( p.z() - zTopCut );
314 
315  if ( (distZBottom < distR) || (distZTop < distR) )
316  {
317  return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
318  }
319  return ( norm *= radius );
320 }
321 
323 //
324 // Calculate distance to shape from outside, along normalised vector
325 // - return kInfinity if no intersection, or intersection distance <= tolerance
326 //
327 
329  const G4ThreeVector& v ) const
330 {
331  G4double distMin = std::min(xSemiAxis,ySemiAxis);
332  const G4double dRmax = 100.*std::min(distMin,zSemiAxis);
333  distMin= kInfinity;
334 
335  // check to see if Z plane is relevant
336  if (p.z() <= zBottomCut+halfCarTolerance)
337  {
338  if (v.z() <= 0.0) { return distMin; }
339  G4double distZ = (zBottomCut - p.z()) / v.z();
340 
341  if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
342  {
343  // early exit since can't intercept curved surface if we reach here
344  if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
345  return distMin= distZ;
346  }
347  }
348  if (p.z() >= zTopCut-halfCarTolerance)
349  {
350  if (v.z() >= 0.0) { return distMin;}
351  G4double distZ = (zTopCut - p.z()) / v.z();
352  if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
353  {
354  // early exit since can't intercept curved surface if we reach here
355  if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
356  return distMin= distZ;
357  }
358  }
359  // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface
360 
361  // now check curved surface intercept
362  G4double A,B,C;
363 
364  A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
365  C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0;
366  B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis)
367  + p.y()*v.y()/(ySemiAxis*ySemiAxis)
368  + p.z()*v.z()/(zSemiAxis*zSemiAxis) );
369 
370  C= B*B - 4.0*A*C;
371  if (C > 0.0)
372  {
373  G4double distR= (-B - std::sqrt(C)) / (2.0*A);
374  G4double intZ = p.z()+distR*v.z();
375  if ( (distR > halfRadTolerance)
376  && (intZ >= zBottomCut-halfRadTolerance)
377  && (intZ <= zTopCut+halfRadTolerance) )
378  {
379  distMin = distR;
380  }
381  else if( (distR >- halfRadTolerance)
382  && (intZ >= zBottomCut-halfRadTolerance)
383  && (intZ <= zTopCut+halfRadTolerance) )
384  {
385  // p is on the curved surface, DistanceToIn returns 0 or kInfinity:
386  // DistanceToIn returns 0, if second root is positive (means going inside)
387  // If second root is negative, DistanceToIn returns kInfinity (outside)
388  //
389  distR = (-B + std::sqrt(C) ) / (2.0*A);
390  if(distR>0.) { distMin=0.; }
391  }
392  else
393  {
394  distR= (-B + std::sqrt(C)) / (2.0*A);
395  intZ = p.z()+distR*v.z();
396  if ( (distR > halfRadTolerance)
397  && (intZ >= zBottomCut-halfRadTolerance)
398  && (intZ <= zTopCut+halfRadTolerance) )
399  {
401  if (norm.dot(v)<0.) { distMin = distR; }
402  }
403  }
404  if ( (distMin!=kInfinity) && (distMin>dRmax) )
405  { // Avoid rounding errors due to precision issues on
406  // 64 bits systems. Split long distances and recompute
407  G4double fTerm = distMin-std::fmod(distMin,dRmax);
408  distMin = fTerm + DistanceToIn(p+fTerm*v,v);
409  }
410  }
411 
412  if (std::fabs(distMin)<halfRadTolerance) { distMin=0.; }
413  return distMin;
414 }
415 
417 //
418 // Calculate distance (<= actual) to closest surface of shape from outside
419 // - Return 0 if point inside
420 
422 {
423  G4double distR, distZ;
424 
425  // normal vector: parallel to normal, magnitude 1/(characteristic radius)
426  //
427  G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
428  p.y()/(ySemiAxis*ySemiAxis),
429  p.z()/(zSemiAxis*zSemiAxis));
430  G4double radius= 1.0/norm.mag();
431 
432  // approximate distance to curved surface ( <= actual distance )
433  //
434  distR= (p*norm - 1.0) * radius / 2.0;
435 
436  // Distance to z-cut plane
437  //
438  distZ= zBottomCut - p.z();
439  if (distZ < 0.0)
440  {
441  distZ = p.z() - zTopCut;
442  }
443 
444  // Distance to closest surface from outside
445  //
446  if (distZ < 0.0)
447  {
448  return (distR < 0.0) ? 0.0 : distR;
449  }
450  else if (distR < 0.0)
451  {
452  return distZ;
453  }
454  else
455  {
456  return (distZ < distR) ? distZ : distR;
457  }
458 }
459 
461 //
462 // Calculate distance to surface of shape from `inside', allowing for tolerance
463 
465  const G4ThreeVector& v,
466  const G4bool calcNorm,
467  G4bool *validNorm,
468  G4ThreeVector *n ) const
469 {
470  G4double distMin;
471  enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
472 
473  distMin= kInfinity;
474  surface= kNoSurf;
475 
476  // check to see if Z plane is relevant
477  //
478  if (v.z() < 0.0)
479  {
480  G4double distZ = (zBottomCut - p.z()) / v.z();
481  if (distZ < 0.0)
482  {
483  distZ= 0.0;
484  if (!calcNorm) {return 0.0;}
485  }
486  distMin= distZ;
487  surface= kPlaneSurf;
488  }
489  if (v.z() > 0.0)
490  {
491  G4double distZ = (zTopCut - p.z()) / v.z();
492  if (distZ < 0.0)
493  {
494  distZ= 0.0;
495  if (!calcNorm) {return 0.0;}
496  }
497  distMin= distZ;
498  surface= kPlaneSurf;
499  }
500 
501  // normal vector: parallel to normal, magnitude 1/(characteristic radius)
502  //
503  G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis),
504  p.y()/(ySemiAxis*ySemiAxis),
505  p.z()/(zSemiAxis*zSemiAxis));
506 
507  // now check curved surface intercept
508  //
509  G4double A,B,C;
510 
511  A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
512  C= (p * nearnorm) - 1.0;
513  B= 2.0 * (v * nearnorm);
514 
515  C= B*B - 4.0*A*C;
516  if (C > 0.0)
517  {
518  G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
519  if (distR < 0.0)
520  {
521  distR= 0.0;
522  if (!calcNorm) {return 0.0;}
523  }
524  if (distR < distMin)
525  {
526  distMin= distR;
527  surface= kCurvedSurf;
528  }
529  }
530 
531  // set normal if requested
532  //
533  if (calcNorm)
534  {
535  if (surface == kNoSurf)
536  {
537  *validNorm = false;
538  }
539  else
540  {
541  *validNorm = true;
542  switch (surface)
543  {
544  case kPlaneSurf:
545  *n= G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
546  break;
547  case kCurvedSurf:
548  {
549  G4ThreeVector pexit= p + distMin*v;
550  G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis),
551  pexit.y()/(ySemiAxis*ySemiAxis),
552  pexit.z()/(zSemiAxis*zSemiAxis));
553  truenorm *= 1.0/truenorm.mag();
554  *n= truenorm;
555  } break;
556  default: // Should never reach this case ...
557  DumpInfo();
558  std::ostringstream message;
559  G4int oldprc = message.precision(16);
560  message << "Undefined side for valid surface normal to solid."
561  << G4endl
562  << "Position:" << G4endl
563  << " p.x() = " << p.x()/mm << " mm" << G4endl
564  << " p.y() = " << p.y()/mm << " mm" << G4endl
565  << " p.z() = " << p.z()/mm << " mm" << G4endl
566  << "Direction:" << G4endl << G4endl
567  << " v.x() = " << v.x() << G4endl
568  << " v.y() = " << v.y() << G4endl
569  << " v.z() = " << v.z() << G4endl
570  << "Proposed distance :" << G4endl
571  << " distMin = " << distMin/mm << " mm";
572  message.precision(oldprc);
573  G4Exception("G4Ellipsoid::DistanceToOut(p,v,..)",
574  "GeomSolids1002", JustWarning, message);
575  break;
576  }
577  }
578  }
579 
580  return distMin;
581 }
582 
584 //
585 // Calculate distance (<=actual) to closest surface of shape from inside
586 
588 {
589  G4double distR, distZ;
590 
591 #ifdef G4SPECSDEBUG
592  if( Inside(p) == kOutside )
593  {
594  DumpInfo();
595  std::ostringstream message;
596  G4int oldprc = message.precision(16);
597  message << "Point p is outside !?" << G4endl
598  << "Position:" << G4endl
599  << " p.x() = " << p.x()/mm << " mm" << G4endl
600  << " p.y() = " << p.y()/mm << " mm" << G4endl
601  << " p.z() = " << p.z()/mm << " mm";
602  message.precision(oldprc) ;
603  G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
604  JustWarning, message);
605  }
606 #endif
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 
614  // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag())
615  //
616  G4double radius= p.mag();
617  G4double tmp= norm.mag();
618  if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
619 
620  // Approximate distance to curved surface ( <= actual distance )
621  //
622  distR = (1.0 - p*norm) * radius / 2.0;
623 
624  // Distance to z-cut plane
625  //
626  distZ = p.z() - zBottomCut;
627  if (distZ < 0.0) {distZ= zTopCut - p.z();}
628 
629  // Distance to closest surface from inside
630  //
631  if ( (distZ < 0.0) || (distR < 0.0) )
632  {
633  return 0.0;
634  }
635  else
636  {
637  return (distZ < distR) ? distZ : distR;
638  }
639 }
640 
642 //
643 // G4EntityType
644 
646 {
647  return G4String("G4Ellipsoid");
648 }
649 
651 //
652 // Make a clone of the object
653 
655 {
656  return new G4Ellipsoid(*this);
657 }
658 
660 //
661 // Stream object contents to an output stream
662 
663 std::ostream& G4Ellipsoid::StreamInfo( std::ostream& os ) const
664 {
665  G4int oldprc = os.precision(16);
666  os << "-----------------------------------------------------------\n"
667  << " *** Dump for solid - " << GetName() << " ***\n"
668  << " ===================================================\n"
669  << " Solid type: G4Ellipsoid\n"
670  << " Parameters: \n"
671 
672  << " semi-axis x: " << xSemiAxis/mm << " mm \n"
673  << " semi-axis y: " << ySemiAxis/mm << " mm \n"
674  << " semi-axis z: " << zSemiAxis/mm << " mm \n"
675  << " max semi-axis: " << semiAxisMax/mm << " mm \n"
676  << " lower cut plane level z: " << zBottomCut/mm << " mm \n"
677  << " upper cut plane level z: " << zTopCut/mm << " mm \n"
678  << "-----------------------------------------------------------\n";
679  os.precision(oldprc);
680 
681  return os;
682 }
683 
685 //
686 // GetPointOnSurface
687 
689 {
690  G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
691  G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3;
692 
693  max1 = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
694  max1 = max1 > zSemiAxis ? max1 : zSemiAxis;
695  if (max1 == xSemiAxis) { max2 = ySemiAxis; max3 = zSemiAxis; }
696  else if (max1 == ySemiAxis) { max2 = xSemiAxis; max3 = zSemiAxis; }
697  else { max2 = xSemiAxis; max3 = ySemiAxis; }
698 
699  phi = G4RandFlat::shoot(0.,twopi);
700 
701  cosphi = std::cos(phi); sinphi = std::sin(phi);
702  costheta = G4RandFlat::shoot(zBottomCut,zTopCut)/zSemiAxis;
703  sintheta = std::sqrt(1.-sqr(costheta));
704 
705  alpha = 1.-sqr(max2/max1); beta = 1.-sqr(max3/max1);
706 
707  aTop = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis));
708  aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis));
709 
710  // approximation
711  // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf"
712  aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)-
713  1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta)));
714 
715  aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis);
716 
717  if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis )
718  || ( zTopCut == 0 && zBottomCut ==0 ) )
719  {
720  aTop = 0; aBottom = 0;
721  }
722 
723  chose = G4RandFlat::shoot(0.,aTop + aBottom + aCurved);
724 
725  if(chose < aCurved)
726  {
727  xRand = xSemiAxis*sintheta*cosphi;
728  yRand = ySemiAxis*sintheta*sinphi;
729  zRand = zSemiAxis*costheta;
730  return G4ThreeVector (xRand,yRand,zRand);
731  }
732  else if(chose >= aCurved && chose < aCurved + aTop)
733  {
734  xRand = G4RandFlat::shoot(-1.,1.)*xSemiAxis
735  * std::sqrt(1-sqr(zTopCut/zSemiAxis));
736  yRand = G4RandFlat::shoot(-1.,1.)*ySemiAxis
737  * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis));
738  zRand = zTopCut;
739  return G4ThreeVector (xRand,yRand,zRand);
740  }
741  else
742  {
743  xRand = G4RandFlat::shoot(-1.,1.)*xSemiAxis
744  * std::sqrt(1-sqr(zBottomCut/zSemiAxis));
745  yRand = G4RandFlat::shoot(-1.,1.)*ySemiAxis
746  * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis));
747  zRand = zBottomCut;
748  return G4ThreeVector (xRand,yRand,zRand);
749  }
750 }
751 
753 //
754 // Methods for visualisation
755 
757 {
758  scene.AddSolid(*this);
759 }
760 
762 {
763  // Define the sides of the box into which the G4Ellipsoid instance would fit.
764  //
765  return G4VisExtent (-semiAxisMax, semiAxisMax,
766  -semiAxisMax, semiAxisMax,
767  -semiAxisMax, semiAxisMax);
768 }
769 
771 {
772  return new G4PolyhedronEllipsoid(xSemiAxis, ySemiAxis, zSemiAxis,
773  zBottomCut, zTopCut);
774 }
775 
777 {
778  if (!fpPolyhedron ||
782  {
783  G4AutoLock l(&polyhedronMutex);
784  delete fpPolyhedron;
786  fRebuildPolyhedron = false;
787  l.unlock();
788  }
789  return fpPolyhedron;
790 }
void set(double x, double y, double z)
G4Polyhedron * GetPolyhedron() const
Definition: G4Ellipsoid.cc:776
G4String GetName() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4GeometryType GetEntityType() const
Definition: G4Ellipsoid.cc:645
G4bool fRebuildPolyhedron
Definition: G4Ellipsoid.hh:134
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double x() const
double dot(const Hep3Vector &) const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Ellipsoid.cc:202
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Ellipsoid.cc:464
const char * p
Definition: xmltok.h:285
G4double GetZTopCut() const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Ellipsoid.cc:252
G4double GetZBottomCut() const
double C(double temp)
double B(double temperature)
virtual ~G4Ellipsoid()
Definition: G4Ellipsoid.cc:135
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:231
int G4int
Definition: G4Types.hh:78
G4double GetSemiAxisMax(G4int i) const
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
double z() const
void DumpInfo() const
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetZCuts(G4double newzBottomCut, G4double newzTopCut)
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Ellipsoid.cc:756
G4Ellipsoid(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double pzSemiAxis, G4double pzBottomCut=0, G4double pzTopCut=0)
Definition: G4Ellipsoid.cc:73
double A(double temperature)
G4VisExtent GetExtent() const
Definition: G4Ellipsoid.cc:761
G4Polyhedron * fpPolyhedron
Definition: G4Ellipsoid.hh:135
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4Ellipsoid & operator=(const G4Ellipsoid &rhs)
Definition: G4Ellipsoid.cc:161
G4Polyhedron * CreatePolyhedron() const
Definition: G4Ellipsoid.cc:770
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4VSolid * Clone() const
Definition: G4Ellipsoid.cc:654
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Ellipsoid.cc:328
const G4int n
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Ellipsoid.cc:294
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
G4int G4Mutex
Definition: G4Threading.hh:173
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Ellipsoid.cc:663
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
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
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:111
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
static constexpr double pi
Definition: G4SIunits.hh:75
T sqr(const T &x)
Definition: templates.hh:145
G4ThreeVector GetPointOnSurface() const
Definition: G4Ellipsoid.cc:688
double G4double
Definition: G4Types.hh:76
static const G4double alpha
double mag() const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Ellipsoid.cc:191
static G4GeometryTolerance * GetInstance()
static int max3(int a, int b, int c)