Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4Ellipsoid.cc
이 파일의 문서화 페이지로 가기
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 104316 2017-05-24 13:04:23Z 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: reimplemented CalculateExtent() using
37 // G4BoundingEnvelope, removed CreateRotatedVertices()
38 //
39 // --------------------------------------------------------------------
40 
41 #include "globals.hh"
42 
43 #include "G4Ellipsoid.hh"
44 
45 #include "G4VoxelLimits.hh"
46 #include "G4AffineTransform.hh"
47 #include "G4GeometryTolerance.hh"
48 #include "G4BoundingEnvelope.hh"
49 
50 #include "meshdefs.hh"
51 #include "Randomize.hh"
52 
53 #include "G4VPVParameterisation.hh"
54 
55 #include "G4VGraphicsScene.hh"
56 #include "G4VisExtent.hh"
57 
58 #include "G4AutoLock.hh"
59 
60 namespace
61 {
62  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
63 }
64 
65 using namespace CLHEP;
66 
68 //
69 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
70 // - note if pDPhi>2PI then reset to 2PI
71 
73  G4double pxSemiAxis,
74  G4double pySemiAxis,
75  G4double pzSemiAxis,
76  G4double pzBottomCut,
77  G4double pzTopCut)
78  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
79  fCubicVolume(0.), fSurfaceArea(0.), zBottomCut(0.), zTopCut(0.)
80 {
81  // note: for users that want to use the full ellipsoid it is useful
82  // to include a default for the cuts
83 
85 
88 
89  // Check Semi-Axis
90  if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
91  {
92  std::ostringstream message;
93  message << "Invalid semi-axis - " << GetName();
94  G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
95  FatalErrorInArgument, message);
96  }
97  SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis);
98 
99  if ( pzBottomCut == 0 && pzTopCut == 0 )
100  {
101  SetZCuts(-pzSemiAxis, pzSemiAxis);
102  }
103  else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
104  && (pzBottomCut < pzTopCut) )
105  {
106  SetZCuts(pzBottomCut, pzTopCut);
107  }
108  else
109  {
110  std::ostringstream message;
111  message << "Invalid z-coordinate for cutting plane - " << GetName();
112  G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
113  FatalErrorInArgument, message);
114  }
115 }
116 
118 //
119 // Fake default constructor - sets only member data and allocates memory
120 // for usage restricted to object persistency.
121 //
123  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0), kRadTolerance(0.),
124  halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
125  fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
126  semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
127 {
128 }
129 
131 //
132 // Destructor
133 
135 {
136  delete fpPolyhedron; fpPolyhedron = 0;
137 }
138 
140 //
141 // Copy constructor
142 
144  : G4VSolid(rhs),
145  fRebuildPolyhedron(false), fpPolyhedron(0),
146  kRadTolerance(rhs.kRadTolerance),
147  halfCarTolerance(rhs.halfCarTolerance),
148  halfRadTolerance(rhs.halfRadTolerance),
149  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
150  xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
151  zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
152  zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
153 {
154 }
155 
157 //
158 // Assignment operator
159 
161 {
162  // Check assignment to self
163  //
164  if (this == &rhs) { return *this; }
165 
166  // Copy base class data
167  //
168  G4VSolid::operator=(rhs);
169 
170  // Copy data
171  //
176  xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
178  zBottomCut = rhs.zBottomCut; zTopCut = rhs.zTopCut;
179  fRebuildPolyhedron = false;
180  delete fpPolyhedron; fpPolyhedron = 0;
181 
182  return *this;
183 }
184 
186 //
187 // Dispatch to parameterisation for replication mechanism dimension
188 // computation & modification.
189 
191  const G4int n,
192  const G4VPhysicalVolume* pRep)
193 {
194  p->ComputeDimensions(*this,n,pRep);
195 }
196 
198 //
199 // Get bounding box
200 
202 {
203  G4double dx = GetSemiAxisMax(0);
204  G4double dy = GetSemiAxisMax(1);
205  G4double dz = GetSemiAxisMax(2);
206  G4double zmin = std::max(-dz,GetZBottomCut());
207  G4double zmax = std::min( dz,GetZTopCut());
208  pMin.set(-dx,-dy,zmin);
209  pMax.set( dx, dy,zmax);
210 
211  // Check correctness of the bounding box
212  //
213  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
214  {
215  std::ostringstream message;
216  message << "Bad bounding box (min >= max) for solid: "
217  << GetName() << " !"
218  << "\npMin = " << pMin
219  << "\npMax = " << pMax;
220  G4Exception("G4Ellipsoid::BoundingLimits()", "GeomMgt0001",
221  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  BoundingLimits(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))
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)
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  //
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 {
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  //
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  //
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);
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  //
768 }
769 
771 {
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 }
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4VSolid * Clone() const
Definition: G4Ellipsoid.cc:654
CLHEP::Hep3Vector G4ThreeVector
G4double GetSemiAxisMax(G4int i) const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
Definition: geomdefs.hh:42
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Ellipsoid.cc:201
G4bool fRebuildPolyhedron
Definition: G4Ellipsoid.hh:134
G4GeometryType GetEntityType() const
Definition: G4Ellipsoid.cc:645
static constexpr double mm
Definition: G4SIunits.hh:115
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
#define G4endl
Definition: G4ios.hh:61
const char * p
Definition: xmltok.h:285
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4double zSemiAxis
Definition: G4Ellipsoid.hh:144
double z() const
double dot(const Hep3Vector &) const
G4double GetRadialTolerance() const
G4VisExtent GetExtent() const
Definition: G4Ellipsoid.cc:761
Float_t tmp
Double_t beta
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Ellipsoid.cc:663
void SetSemiAxis(G4double x, G4double y, G4double z)
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
G4double zBottomCut
Definition: G4Ellipsoid.hh:144
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Ellipsoid.cc:231
static int max3(int a, int b, int c)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Ellipsoid.cc:294
G4double halfCarTolerance
Definition: G4Ellipsoid.hh:140
G4ThreeVector GetPointOnSurface() const
Definition: G4Ellipsoid.cc:688
virtual ~G4Ellipsoid()
Definition: G4Ellipsoid.cc:134
G4double halfRadTolerance
Definition: G4Ellipsoid.hh:140
G4Polyhedron * GetPolyhedron() const
Definition: G4Ellipsoid.cc:776
static const G4double alpha
G4double kCarTolerance
Definition: G4VSolid.hh:307
static constexpr double twopi
Definition: G4SIunits.hh:76
double A(double temperature)
G4double fCubicVolume
Definition: G4Ellipsoid.hh:142
G4double ySemiAxis
Definition: G4Ellipsoid.hh:144
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Ellipsoid.cc:328
G4double xSemiAxis
Definition: G4Ellipsoid.hh:144
Double_t radius
G4double semiAxisMax
Definition: G4Ellipsoid.hh:144
G4double GetZBottomCut() const
void SetZCuts(G4double newzBottomCut, G4double newzTopCut)
ThreeVector shoot(const G4int Ap, const G4int Af)
static G4int GetNumberOfRotationSteps()
G4double fSurfaceArea
Definition: G4Ellipsoid.hh:143
G4Ellipsoid(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double pzSemiAxis, G4double pzBottomCut=0, G4double pzTopCut=0)
Definition: G4Ellipsoid.cc:72
EInside Inside(const G4ThreeVector &p) const
Definition: G4Ellipsoid.cc:252
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
double mag() const
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
int G4int
Definition: G4Types.hh:78
EInside
Definition: geomdefs.hh:58
ifstream in
Definition: comparison.C:7
G4Polyhedron * fpPolyhedron
Definition: G4Ellipsoid.hh:135
EAxis
Definition: geomdefs.hh:54
G4String GetName() const
void DumpInfo() const
Float_t norm
double C(double temp)
G4Polyhedron * CreatePolyhedron() const
Definition: G4Ellipsoid.cc:770
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Ellipsoid.cc:756
double x() const
G4double GetZTopCut() const
G4double kRadTolerance
Definition: G4Ellipsoid.hh:139
static G4GeometryTolerance * GetInstance()
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Ellipsoid.cc:190
Char_t n[5]
T sqr(const T &x)
Definition: templates.hh:145
static constexpr double pi
Definition: G4SIunits.hh:75
double y() const
G4Ellipsoid & operator=(const G4Ellipsoid &rhs)
Definition: G4Ellipsoid.cc:160
G4double zTopCut
Definition: G4Ellipsoid.hh:144
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Ellipsoid.cc:464
double B(double temperature)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
std::mutex G4Mutex
Definition: G4Threading.hh:84