ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EllipticalCone.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EllipticalCone.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 // Implementation of G4EllipticalCone class
27 //
28 // This code implements an Elliptical Cone given explicitly by the
29 // equation:
30 // x^2/a^2 + y^2/b^2 = (z-h)^2
31 // and specified by the parameters (a,b,h) and a cut parallel to the
32 // xy plane above z = 0.
33 //
34 // Author: Dionysios Anninos
35 // Revised: Evgueni Tcherniaev
36 // --------------------------------------------------------------------
37 
38 #if !(defined(G4GEOM_USE_UELLIPTICALCONE) && defined(G4GEOM_USE_SYS_USOLIDS))
39 
40 #include "globals.hh"
41 
42 #include "G4EllipticalCone.hh"
43 
44 #include "G4RandomTools.hh"
45 #include "G4GeomTools.hh"
46 #include "G4ClippablePolygon.hh"
47 #include "G4VoxelLimits.hh"
48 #include "G4AffineTransform.hh"
49 #include "G4BoundingEnvelope.hh"
50 #include "G4GeometryTolerance.hh"
51 
52 #include "meshdefs.hh"
53 
54 #include "Randomize.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
71 
73  G4double pxSemiAxis,
74  G4double pySemiAxis,
75  G4double pzMax,
76  G4double pzTopCut)
77  : G4VSolid(pName), zTopCut(0.)
78 {
80 
81  // Check Semi-Axis & Z-cut
82  //
83  if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) )
84  {
85  std::ostringstream message;
86  message << "Invalid semi-axis or height for solid: " << GetName()
87  << "\n X semi-axis, Y semi-axis, height = "
88  << pxSemiAxis << ", " << pySemiAxis << ", " << pzMax;
89  G4Exception("G4EllipticalCone::G4EllipticalCone()", "GeomSolids0002",
90  FatalErrorInArgument, message);
91  }
92 
93  if ( pzTopCut <= 0 )
94  {
95  std::ostringstream message;
96  message << "Invalid z-coordinate for cutting plane for solid: " << GetName()
97  << "\n Z top cut = " << pzTopCut;
98  G4Exception("G4EllipticalCone::G4EllipticalCone()", "GeomSolids0002",
99  FatalErrorInArgument, message);
100  }
101 
102  SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax );
103  SetZCut(pzTopCut);
104 }
105 
107 //
108 // Fake default constructor - sets only member data and allocates memory
109 // for usage restricted to object persistency.
110 
112  : G4VSolid(a), halfCarTol(0.),
113  xSemiAxis(0.), ySemiAxis(0.), zheight(0.), zTopCut(0.),
114  cosAxisMin(0.), invXX(0.), invYY(0.)
115 {
116 }
117 
119 //
120 // Destructor
121 
123 {
124  delete fpPolyhedron; fpPolyhedron = nullptr;
125 }
126 
128 //
129 // Copy constructor
130 
132  : G4VSolid(rhs), halfCarTol(rhs.halfCarTol),
133  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
134  xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
135  zheight(rhs.zheight), zTopCut(rhs.zTopCut),
136  cosAxisMin(rhs.cosAxisMin), invXX(rhs.invXX), invYY(rhs.invYY)
137 {
138 }
139 
141 //
142 // Assignment operator
143 
145 {
146  // Check assignment to self
147  //
148  if (this == &rhs) { return *this; }
149 
150  // Copy base class data
151  //
152  G4VSolid::operator=(rhs);
153 
154  // Copy data
155  //
156  halfCarTol = rhs.halfCarTol;
158  xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
159  zheight = rhs.zheight; zTopCut = rhs.zTopCut;
160  cosAxisMin = rhs.cosAxisMin; invXX = rhs.invXX; invYY = rhs.invYY;
161 
162  fRebuildPolyhedron = false;
163  delete fpPolyhedron; fpPolyhedron = nullptr;
164 
165  return *this;
166 }
167 
169 //
170 // Get bounding box
171 
173  G4ThreeVector& pMax) const
174 {
175  G4double zcut = GetZTopCut();
176  G4double height = GetZMax();
177  G4double xmax = GetSemiAxisX()*(height+zcut);
178  G4double ymax = GetSemiAxisY()*(height+zcut);
179  pMin.set(-xmax,-ymax,-zcut);
180  pMax.set( xmax, ymax, zcut);
181 
182  // Check correctness of the bounding box
183  //
184  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
185  {
186  std::ostringstream message;
187  message << "Bad bounding box (min >= max) for solid: "
188  << GetName() << " !"
189  << "\npMin = " << pMin
190  << "\npMax = " << pMax;
191  G4Exception("G4EllipticalCone::BoundingLimits()", "GeomMgt0001",
192  JustWarning, message);
193  DumpInfo();
194  }
195 }
196 
198 //
199 // Calculate extent under transform and specified limit
200 
201 G4bool
203  const G4VoxelLimits& pVoxelLimit,
204  const G4AffineTransform& pTransform,
205  G4double& pMin, G4double& pMax) const
206 {
207  G4ThreeVector bmin,bmax;
208  G4bool exist;
209 
210  // Check bounding box (bbox)
211  //
212  BoundingLimits(bmin,bmax);
213  G4BoundingEnvelope bbox(bmin,bmax);
214 #ifdef G4BBOX_EXTENT
215  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
216 #endif
217  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
218  {
219  return exist = (pMin < pMax) ? true : false;
220  }
221 
222  // Set bounding envelope (benv) and calculate extent
223  //
224  static const G4int NSTEPS = 48; // number of steps for whole circle
225  static const G4double ang = twopi/NSTEPS;
226  static const G4double sinHalf = std::sin(0.5*ang);
227  static const G4double cosHalf = std::cos(0.5*ang);
228  static const G4double sinStep = 2.*sinHalf*cosHalf;
229  static const G4double cosStep = 1. - 2.*sinHalf*sinHalf;
230  G4double zcut = bmax.z();
231  G4double height = GetZMax();
232  G4double sxmin = GetSemiAxisX()*(height-zcut)/cosHalf;
233  G4double symin = GetSemiAxisY()*(height-zcut)/cosHalf;
234  G4double sxmax = bmax.x()/cosHalf;
235  G4double symax = bmax.y()/cosHalf;
236 
237  G4double sinCur = sinHalf;
238  G4double cosCur = cosHalf;
239  G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
240  for (G4int k=0; k<NSTEPS; ++k)
241  {
242  baseA[k].set(sxmax*cosCur,symax*sinCur,-zcut);
243  baseB[k].set(sxmin*cosCur,symin*sinCur, zcut);
244 
245  G4double sinTmp = sinCur;
246  sinCur = sinCur*cosStep + cosCur*sinStep;
247  cosCur = cosCur*cosStep - sinTmp*sinStep;
248  }
249 
250  std::vector<const G4ThreeVectorList *> polygons(2);
251  polygons[0] = &baseA;
252  polygons[1] = &baseB;
253  G4BoundingEnvelope benv(bmin,bmax,polygons);
254  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
255  return exist;
256 }
257 
259 //
260 // Determine where is point: inside, outside or on surface
261 
263 {
264  G4double hp = std::sqrt(p.x()*p.x()*invXX + p.y()*p.y()*invYY) + p.z();
265  G4double ds = (hp - zheight)*cosAxisMin;
266  G4double dz = std::abs(p.z()) - zTopCut;
267  G4double dist = std::max(ds,dz);
268 
269  if (dist > halfCarTol) return kOutside;
270  return (dist > -halfCarTol) ? kSurface : kInside;
271 }
272 
274 //
275 // Return unit normal at surface closest to p
276 
278 {
279  G4ThreeVector norm(0,0,0);
280  G4int nsurf = 0; // number of surfaces where p is placed
281 
282  G4double hp = std::sqrt(p.x()*p.x()*invXX + p.y()*p.y()*invYY) + p.z();
283  G4double ds = (hp - zheight)*cosAxisMin;
284  if (std::abs(ds) <= halfCarTol)
285  {
286  norm = G4ThreeVector(p.x()*invXX, p.y()*invYY, hp - p.z());
287  G4double mag = norm.mag();
288  if (mag == 0) return G4ThreeVector(0,0,1); // apex
289  norm *= (1/mag);
290  ++nsurf;
291  }
292  G4double dz = std::abs(p.z()) - zTopCut;
293  if (std::abs(dz) <= halfCarTol)
294  {
295  norm += G4ThreeVector(0., 0.,(p.z() < 0) ? -1. : 1.);
296  ++nsurf;
297  }
298 
299  if (nsurf == 1) return norm;
300  else if (nsurf > 1) return norm.unit(); // elliptic edge
301  else
302  {
303  // Point is not on the surface
304  //
305 #ifdef G4CSGDEBUG
306  std::ostringstream message;
307  G4int oldprc = message.precision(16);
308  message << "Point p is not on surface (!?) of solid: "
309  << GetName() << G4endl;
310  message << "Position:\n";
311  message << " p.x() = " << p.x()/mm << " mm\n";
312  message << " p.y() = " << p.y()/mm << " mm\n";
313  message << " p.z() = " << p.z()/mm << " mm";
314  G4cout.precision(oldprc);
315  G4Exception("G4EllipticalCone::SurfaceNormal(p)", "GeomSolids1002",
316  JustWarning, message );
317  DumpInfo();
318 #endif
319  return ApproxSurfaceNormal(p);
320  }
321 }
322 
324 //
325 // Find surface nearest to point and return corresponding normal.
326 // The algorithm is similar to the algorithm used in Inside().
327 // This method normally should not be called.
328 
331 {
332  G4double hp = std::sqrt(p.x()*p.x()*invXX + p.y()*p.y()*invYY) + p.z();
333  G4double ds = (hp - zheight)*cosAxisMin;
334  G4double dz = std::abs(p.z()) - zTopCut;
335  if (ds > dz && std::abs(hp - p.z()) > halfCarTol)
336  return G4ThreeVector(p.x()*invXX, p.y()*invYY, hp - p.z()).unit();
337  else
338  return G4ThreeVector(0., 0.,(p.z() < 0) ? -1. : 1.);
339 }
340 
342 //
343 // Calculate distance to shape from outside, along normalised vector
344 // return kInfinity if no intersection, or intersection distance <= tolerance
345 
347  const G4ThreeVector& v ) const
348 {
349  G4double distMin = kInfinity;
350 
351  // code from EllipticalTube
352 
353  G4double sigz = p.z()+zTopCut;
354 
355  //
356  // Check z = -dz planer surface
357  //
358 
359  if (sigz < halfCarTol)
360  {
361  //
362  // We are "behind" the shape in z, and so can
363  // potentially hit the rear face. Correct direction?
364  //
365  if (v.z() <= 0)
366  {
367  //
368  // As long as we are far enough away, we know we
369  // can't intersect
370  //
371  if (sigz < 0) return kInfinity;
372 
373  //
374  // Otherwise, we don't intersect unless we are
375  // on the surface of the ellipse
376  //
377 
378  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
379  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight + zTopCut ) )
380  return kInfinity;
381 
382  }
383  else
384  {
385  //
386  // How far?
387  //
388  G4double q = -sigz/v.z();
389 
390  //
391  // Where does that place us?
392  //
393  G4double xi = p.x() + q*v.x(),
394  yi = p.y() + q*v.y();
395 
396  //
397  // Is this on the surface (within ellipse)?
398  //
399  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) )
400  {
401  //
402  // Yup. Return q, unless we are on the surface
403  //
404  return (sigz < -halfCarTol) ? q : 0;
405  }
406  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
407  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
408  {
409  //
410  // Else, if we are traveling outwards, we know
411  // we must miss
412  //
413  // return kInfinity;
414  }
415  }
416  }
417 
418  //
419  // Check z = +dz planer surface
420  //
421  sigz = p.z() - zTopCut;
422 
423  if (sigz > -halfCarTol)
424  {
425  if (v.z() >= 0)
426  {
427 
428  if (sigz > 0) return kInfinity;
429 
430  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
431  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight-zTopCut ) )
432  return kInfinity;
433 
434  }
435  else {
436  G4double q = -sigz/v.z();
437 
438  G4double xi = p.x() + q*v.x(),
439  yi = p.y() + q*v.y();
440 
441  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) )
442  {
443  return (sigz > -halfCarTol) ? q : 0;
444  }
445  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
446  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
447  {
448  // return kInfinity;
449  }
450  }
451  }
452 
453 
454 #if 0
455 
456  // check to see if Z plane is relevant
457  //
458  if (p.z() < -zTopCut - halfCarTol)
459  {
460  if (v.z() <= 0.0)
461  return distMin;
462 
463  G4double lambda = (-zTopCut - p.z())/v.z();
464 
465  if ( sqr((lambda*v.x()+p.x())/xSemiAxis) +
466  sqr((lambda*v.y()+p.y())/ySemiAxis) <=
467  sqr(zTopCut + zheight + halfCarTol) )
468  {
469  return distMin = std::fabs(lambda);
470  }
471  }
472 
473  if (p.z() > zTopCut + halfCarTol)
474  {
475  if (v.z() >= 0.0)
476  { return distMin; }
477 
478  G4double lambda = (zTopCut - p.z()) / v.z();
479 
480  if ( sqr((lambda*v.x() + p.x())/xSemiAxis) +
481  sqr((lambda*v.y() + p.y())/ySemiAxis) <=
483  {
484  return distMin = std::fabs(lambda);
485  }
486  }
487 
488  if (p.z() > zTopCut - halfCarTol
489  && p.z() < zTopCut + halfCarTol )
490  {
491  if (v.z() > 0.)
492  { return kInfinity; }
493 
494  return distMin = 0.;
495  }
496 
497  if (p.z() < -zTopCut + halfCarTol
498  && p.z() > -zTopCut - halfCarTol)
499  {
500  if (v.z() < 0.)
501  { return distMin = kInfinity; }
502 
503  return distMin = 0.;
504  }
505 
506 #endif
507 
508  // if we are here then it either intersects or grazes the curved surface
509  // or it does not intersect at all
510  //
511  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
512  G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) +
513  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
514  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) -
515  sqr(zheight - p.z());
516 
517  G4double discr = B*B - 4.*A*C;
518 
519  // if the discriminant is negative it never hits the curved object
520  //
521  if ( discr < -halfCarTol )
522  { return distMin; }
523 
524  // case below is when it hits or grazes the surface
525  //
526  if ( (discr >= -halfCarTol ) && (discr < halfCarTol ) )
527  {
528  return distMin = std::fabs(-B/(2.*A));
529  }
530 
531  G4double plus = (-B+std::sqrt(discr))/(2.*A);
532  G4double minus = (-B-std::sqrt(discr))/(2.*A);
533 
534  // Special case::Point on Surface, Check norm.dot(v)
535 
536  if ( ( std::fabs(plus) < halfCarTol )||( std::fabs(minus) < halfCarTol ) )
537  {
538  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
539  p.y()/(ySemiAxis*ySemiAxis),
540  -( p.z() - zheight ));
541  if ( truenorm*v >= 0) // going outside the solid from surface
542  {
543  return kInfinity;
544  }
545  else
546  {
547  return 0;
548  }
549  }
550 
551  // G4double lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
552  G4double lambda = 0;
553 
554  if ( minus > halfCarTol && minus < distMin )
555  {
556  lambda = minus ;
557  // check normal vector n * v < 0
558  G4ThreeVector pin = p + lambda*v;
559  if(std::fabs(pin.z())< zTopCut + halfCarTol)
560  {
561  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
562  pin.y()/(ySemiAxis*ySemiAxis),
563  - ( pin.z() - zheight ));
564  if ( truenorm*v < 0)
565  { // yes, going inside the solid
566  distMin = lambda;
567  }
568  }
569  }
570  if ( plus > halfCarTol && plus < distMin )
571  {
572  lambda = plus ;
573  // check normal vector n * v < 0
574  G4ThreeVector pin = p + lambda*v;
575  if(std::fabs(pin.z()) < zTopCut + halfCarTol)
576  {
577  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
578  pin.y()/(ySemiAxis*ySemiAxis),
579  - ( pin.z() - zheight ) );
580  if ( truenorm*v < 0)
581  { // yes, going inside the solid
582  distMin = lambda;
583  }
584  }
585  }
586  if (distMin < halfCarTol) distMin=0.;
587  return distMin ;
588 }
589 
591 //
592 // Calculate distance (<= actual) to closest surface of shape from outside
593 // Return 0 if point inside
594 
596 {
597  G4double hp = std::sqrt(p.x()*p.x()*invXX + p.y()*p.y()*invYY) + p.z();
598  G4double ds = (hp - zheight)*cosAxisMin;
599  G4double dz = std::abs(p.z()) - zTopCut;
600  G4double dist = std::max(ds,dz);
601  return (dist > 0) ? dist : 0.;
602 }
603 
605 //
606 // Calculate distance to surface of shape from `inside',
607 // allowing for tolerance
608 
610  const G4ThreeVector& v,
611  const G4bool calcNorm,
612  G4bool* validNorm,
613  G4ThreeVector* n ) const
614 {
615  G4double distMin, lambda;
616  enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
617 
618  distMin = kInfinity;
619  surface = kNoSurf;
620 
621  if (v.z() < 0.0)
622  {
623  lambda = (-p.z() - zTopCut)/v.z();
624 
625  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) +
626  sqr((p.y() + lambda*v.y())/ySemiAxis)) <
628  {
629  distMin = std::fabs(lambda);
630 
631  if (!calcNorm) { return distMin; }
632  }
633  distMin = std::fabs(lambda);
634  surface = kPlaneSurf;
635  }
636 
637  if (v.z() > 0.0)
638  {
639  lambda = (zTopCut - p.z()) / v.z();
640 
641  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis)
642  + sqr((p.y() + lambda*v.y())/ySemiAxis) )
643  < (sqr(zheight - zTopCut + halfCarTol)) )
644  {
645  distMin = std::fabs(lambda);
646  if (!calcNorm) { return distMin; }
647  }
648  distMin = std::fabs(lambda);
649  surface = kPlaneSurf;
650  }
651 
652  // if we are here then it either intersects or grazes the
653  // curved surface...
654  //
655  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
656  G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) +
657  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
658  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
659  - sqr(zheight - p.z());
660 
661  G4double discr = B*B - 4.*A*C;
662 
663  if ( discr >= - halfCarTol && discr < halfCarTol )
664  {
665  if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); }
666  }
667 
668  else if ( discr > halfCarTol )
669  {
670  G4double plus = (-B+std::sqrt(discr))/(2.*A);
671  G4double minus = (-B-std::sqrt(discr))/(2.*A);
672 
673  if ( plus > halfCarTol && minus > halfCarTol )
674  {
675  // take the shorter distance
676  //
677  lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
678  }
679  else
680  {
681  // at least one solution is close to zero or negative
682  // so, take small positive solution or zero
683  //
684  lambda = plus > -halfCarTol ? plus : 0;
685  }
686 
687  if ( std::fabs(lambda) < distMin )
688  {
689  if( std::fabs(lambda) > halfCarTol)
690  {
691  distMin = std::fabs(lambda);
692  surface = kCurvedSurf;
693  }
694  else // Point is On the Surface, Check Normal
695  {
696  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
697  p.y()/(ySemiAxis*ySemiAxis),
698  -( p.z() - zheight ));
699  if( truenorm.dot(v) > 0 )
700  {
701  distMin = 0.0;
702  surface = kCurvedSurf;
703  }
704  }
705  }
706  }
707 
708  // set normal if requested
709  //
710  if (calcNorm)
711  {
712  if (surface == kNoSurf)
713  {
714  *validNorm = false;
715  }
716  else
717  {
718  *validNorm = true;
719  switch (surface)
720  {
721  case kPlaneSurf:
722  {
723  *n = G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
724  }
725  break;
726 
727  case kCurvedSurf:
728  {
729  G4ThreeVector pexit = p + distMin*v;
730  G4ThreeVector truenorm( pexit.x()/(xSemiAxis*xSemiAxis),
731  pexit.y()/(ySemiAxis*ySemiAxis),
732  -( pexit.z() - zheight ) );
733  truenorm /= truenorm.mag();
734  *n= truenorm;
735  }
736  break;
737 
738  default: // Should never reach this case ...
739  DumpInfo();
740  std::ostringstream message;
741  G4int oldprc = message.precision(16);
742  message << "Undefined side for valid surface normal to solid."
743  << G4endl
744  << "Position:" << G4endl
745  << " p.x() = " << p.x()/mm << " mm" << G4endl
746  << " p.y() = " << p.y()/mm << " mm" << G4endl
747  << " p.z() = " << p.z()/mm << " mm" << G4endl
748  << "Direction:" << G4endl
749  << " v.x() = " << v.x() << G4endl
750  << " v.y() = " << v.y() << G4endl
751  << " v.z() = " << v.z() << G4endl
752  << "Proposed distance :" << G4endl
753  << " distMin = " << distMin/mm << " mm";
754  message.precision(oldprc);
755  G4Exception("G4EllipticalCone::DistanceToOut(p,v,..)",
756  "GeomSolids1002", JustWarning, message);
757  break;
758  }
759  }
760  }
761 
762  if (distMin < halfCarTol) { distMin=0; }
763 
764  return distMin;
765 }
766 
768 //
769 // Calculate distance (<=actual) to closest surface of shape from inside
770 
772 {
773 #ifdef G4SPECSDEBUG
774  if( Inside(p) == kOutside )
775  {
776  std::ostringstream message;
777  G4int oldprc = message.precision(16);
778  message << "Point p is outside (!?) of solid: " << GetName() << "\n"
779  << "Position:\n"
780  << " p.x() = " << p.x()/mm << " mm\n"
781  << " p.y() = " << p.y()/mm << " mm\n"
782  << " p.z() = " << p.z()/mm << " mm";
783  message.precision(oldprc) ;
784  G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
785  JustWarning, message);
786  DumpInfo();
787  }
788 #endif
789  G4double hp = std::sqrt(p.x()*p.x()*invXX + p.y()*p.y()*invYY) + p.z();
790  G4double ds = (zheight - hp)*cosAxisMin;
791  G4double dz = zTopCut - std::abs(p.z());
792  G4double dist = std::min(ds,dz);
793  return (dist > 0) ? dist : 0.;
794 }
795 
797 //
798 // GetEntityType
799 
801 {
802  return G4String("G4EllipticalCone");
803 }
804 
806 //
807 // Make a clone of the object
808 
810 {
811  return new G4EllipticalCone(*this);
812 }
813 
815 //
816 // Stream object contents to an output stream
817 
818 std::ostream& G4EllipticalCone::StreamInfo( std::ostream& os ) const
819 {
820  G4int oldprc = os.precision(16);
821  os << "-----------------------------------------------------------\n"
822  << " *** Dump for solid - " << GetName() << " ***\n"
823  << " ===================================================\n"
824  << " Solid type: G4EllipticalCone\n"
825  << " Parameters: \n"
826 
827  << " semi-axis x: " << xSemiAxis/mm << " mm \n"
828  << " semi-axis y: " << ySemiAxis/mm << " mm \n"
829  << " height z: " << zheight/mm << " mm \n"
830  << " half length in z: " << zTopCut/mm << " mm \n"
831  << "-----------------------------------------------------------\n";
832  os.precision(oldprc);
833 
834  return os;
835 }
836 
838 //
839 // Return random point on the surface of the solid
840 
842 {
843  G4double x0 = xSemiAxis*zheight; // x semi axis at z=0
844  G4double y0 = ySemiAxis*zheight; // y semi axis at z=0
846  G4double kmin = (zTopCut >= zheight ) ? 0. : (zheight - zTopCut)/zheight;
847  G4double kmax = (zTopCut >= zheight ) ? 2. : (zheight + zTopCut)/zheight;
848 
849  // Set areas (base at -Z, side surface, base at +Z)
850  //
851  G4double szmin = pi*x0*y0*kmax*kmax;
852  G4double szmax = pi*x0*y0*kmin*kmin;
853  G4double sside = s0*(kmax*kmax - kmin*kmin);
854  G4double ssurf[3] = { szmin, sside, szmax };
855  for (auto i=1; i<3; ++i) { ssurf[i] += ssurf[i-1]; }
856 
857  // Select surface
858  //
859  G4double select = ssurf[2]*G4UniformRand();
860  G4int k = 2;
861  if (select <= ssurf[1]) k = 1;
862  if (select <= ssurf[0]) k = 0;
863 
864  // Pick random point on selected surface
865  //
867  switch(k)
868  {
869  case 0: // base at -Z, uniform distribution, rejection sampling
870  {
871  G4double zh = zheight + zTopCut;
873  p.set(rho.x(),rho.y(),-zTopCut);
874  break;
875  }
876  case 1: // side surface, uniform distribution, rejection sampling
877  {
878  G4double zh = G4RandomRadiusInRing(zheight-zTopCut, zheight+zTopCut);
879  G4double a = x0;
880  G4double b = y0;
881 
882  G4double hh = zheight*zheight;
883  G4double aa = a*a;
884  G4double bb = b*b;
885  G4double R = std::max(a,b);
886  G4double mu_max = R*std::sqrt(hh + R*R);
887 
888  G4double x,y;
889  for (auto i=0; i<1000; ++i)
890  {
892  x = std::cos(phi);
893  y = std::sin(phi);
894  G4double xx = x*x;
895  G4double yy = y*y;
896  G4double E = hh + aa*xx + bb*yy;
897  G4double F = (aa-bb)*x*y;
898  G4double G = aa*yy + bb*xx;
899  G4double mu = std::sqrt(E*G - F*F);
900  if (mu_max*G4UniformRand() <= mu) break;
901  }
902  p.set(zh*xSemiAxis*x,zh*ySemiAxis*y,zheight-zh);
903  break;
904  }
905  case 2: // base at +Z, uniform distribution, rejection sampling
906  {
907  G4double zh = zheight - zTopCut;
909  p.set(rho.x(),rho.y(),zTopCut);
910  break;
911  }
912  }
913  return p;
914 }
915 
917 //
918 // Get cubic volume
919 
921 {
922  if (fCubicVolume == 0.0)
923  {
924  G4double x0 = xSemiAxis*zheight; // x semi axis at z=0
925  G4double y0 = ySemiAxis*zheight; // y semi axis at z=0
926  G4double v0 = CLHEP::pi*x0*y0*zheight/3.;
927  G4double kmin = (zTopCut >= zheight ) ? 0. : (zheight - zTopCut)/zheight;
928  G4double kmax = (zTopCut >= zheight ) ? 2. : (zheight + zTopCut)/zheight;
929  fCubicVolume = (kmax - kmin)*(kmax*kmax + kmax*kmin + kmin*kmin)*v0;
930  }
931  return fCubicVolume;
932 }
933 
935 //
936 // Get surface area
937 
939 {
940  if (fSurfaceArea == 0.0)
941  {
942  G4double x0 = xSemiAxis*zheight; // x semi axis at z=0
943  G4double y0 = ySemiAxis*zheight; // y semi axis at z=0
945  G4double kmin = (zTopCut >= zheight ) ? 0. : (zheight - zTopCut)/zheight;
946  G4double kmax = (zTopCut >= zheight ) ? 2. : (zheight + zTopCut)/zheight;
947  fSurfaceArea = (kmax - kmin)*(kmax + kmin)*s0
948  + CLHEP::pi*x0*y0*(kmin*kmin + kmax*kmax);
949  }
950  return fSurfaceArea;
951 }
952 
954 //
955 // Methods for visualisation
956 
958 {
959  scene.AddSolid(*this);
960 }
961 
963 {
964  // Define the sides of the box into which the solid instance would fit.
965  //
966  G4ThreeVector pmin,pmax;
967  BoundingLimits(pmin,pmax);
968  return G4VisExtent(pmin.x(),pmax.x(),
969  pmin.y(),pmax.y(),
970  pmin.z(),pmax.z());
971 }
972 
974 {
976 }
977 
979 {
980  if ( (fpPolyhedron == nullptr)
984  {
985  G4AutoLock l(&polyhedronMutex);
986  delete fpPolyhedron;
988  fRebuildPolyhedron = false;
989  l.unlock();
990  }
991  return fpPolyhedron;
992 }
993 
994 #endif // !defined(G4GEOM_USE_UELLIPTICALCONE) || !defined(G4GEOM_USE_SYS_USOLIDS)