Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4EllipticalCone.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: G4EllipticalCone.cc 101118 2016-11-07 09:10:59Z gcosmo $
27 //
28 // Implementation of G4EllipticalCone class
29 //
30 // This code implements an Elliptical Cone given explicitly by the
31 // equation:
32 // x^2/a^2 + y^2/b^2 = (z-h)^2
33 // and specified by the parameters (a,b,h) and a cut parallel to the
34 // xy plane above z = 0.
35 //
36 // Author: Dionysios Anninos
37 //
38 // --------------------------------------------------------------------
39 
40 #include "globals.hh"
41 
42 #include "G4EllipticalCone.hh"
43 
44 #include "G4ClippablePolygon.hh"
45 #include "G4VoxelLimits.hh"
46 #include "G4AffineTransform.hh"
47 #include "G4BoundingEnvelope.hh"
48 #include "G4GeometryTolerance.hh"
49 
50 #include "meshdefs.hh"
51 
52 #include "Randomize.hh"
53 
54 #include "G4VGraphicsScene.hh"
55 #include "G4VisExtent.hh"
56 
57 #include "G4AutoLock.hh"
58 
59 namespace
60 {
61  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
62 }
63 
64 using namespace CLHEP;
65 
67 //
68 // Constructor - check parameters
69 //
71  G4double pxSemiAxis,
72  G4double pySemiAxis,
73  G4double pzMax,
74  G4double pzTopCut)
75  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
76  fCubicVolume(0.), fSurfaceArea(0.), zTopCut(0.)
77 {
78 
80 
81  halfRadTol = 0.5*kRadTolerance;
82  halfCarTol = 0.5*kCarTolerance;
83 
84  // Check Semi-Axis & Z-cut
85  //
86  if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) )
87  {
88  std::ostringstream message;
89  message << "Invalid semi-axis or height - " << GetName();
90  G4Exception("G4EllipticalCone::G4EllipticalCone()", "GeomSolids0002",
91  FatalErrorInArgument, message);
92  }
93  if ( pzTopCut <= 0 )
94  {
95  std::ostringstream message;
96  message << "Invalid z-coordinate for cutting plane - " << GetName();
97  G4Exception("G4EllipticalCone::G4EllipticalCone()", "InvalidSetup",
98  FatalErrorInArgument, message);
99  }
100 
101  SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax );
102  SetZCut(pzTopCut);
103 }
104 
106 //
107 // Fake default constructor - sets only member data and allocates memory
108 // for usage restricted to object persistency.
109 //
111  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
112  kRadTolerance(0.), halfRadTol(0.), halfCarTol(0.), fCubicVolume(0.),
113  fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zheight(0.),
114  semiAxisMax(0.), zTopCut(0.)
115 {
116 }
117 
119 //
120 // Destructor
121 //
123 {
124  delete fpPolyhedron; fpPolyhedron = 0;
125 }
126 
128 //
129 // Copy constructor
130 //
132  : G4VSolid(rhs),
133  fRebuildPolyhedron(false), fpPolyhedron(0),
134  kRadTolerance(rhs.kRadTolerance),
135  halfRadTol(rhs.halfRadTol), halfCarTol(rhs.halfCarTol),
136  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
137  xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis), zheight(rhs.zheight),
138  semiAxisMax(rhs.semiAxisMax), zTopCut(rhs.zTopCut)
139 {
140 }
141 
143 //
144 // Assignment operator
145 //
147 {
148  // Check assignment to self
149  //
150  if (this == &rhs) { return *this; }
151 
152  // Copy base class data
153  //
154  G4VSolid::operator=(rhs);
155 
156  // Copy data
157  //
158  kRadTolerance = rhs.kRadTolerance;
159  halfRadTol = rhs.halfRadTol; halfCarTol = rhs.halfCarTol;
160  fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
161  xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
162  zheight = rhs.zheight; semiAxisMax = rhs.semiAxisMax; zTopCut = rhs.zTopCut;
163  fRebuildPolyhedron = false;
164  delete fpPolyhedron; fpPolyhedron = 0;
165 
166  return *this;
167 }
168 
170 //
171 // Get bounding box
172 
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::Extent()", "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  Extent(bmin,bmax);
213  G4BoundingEnvelope bbox(bmin,bmax);
214 #ifdef G4BBOX_EXTENT
215  if (true) 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 // Return whether point inside/outside/on surface
261 // Split into radius, phi, theta checks
262 // Each check modifies `in', or returns as approprate
263 //
265 {
266  G4double rad2oo, // outside surface outer tolerance
267  rad2oi; // outside surface inner tolerance
268 
269  EInside in;
270 
271  // check this side of z cut first, because that's fast
272  //
273 
274  if ( (p.z() < -zTopCut - halfCarTol)
275  || (p.z() > zTopCut + halfCarTol ) )
276  {
277  return in = kOutside;
278  }
279 
280  rad2oo= sqr(p.x()/( xSemiAxis + halfRadTol ))
281  + sqr(p.y()/( ySemiAxis + halfRadTol ));
282 
283  if ( rad2oo > sqr( zheight-p.z() ) )
284  {
285  return in = kOutside;
286  }
287 
288  // rad2oi= sqr( p.x()*(1.0 + 0.5*kRadTolerance/(xSemiAxis*xSemiAxis)) )
289  // + sqr( p.y()*(1.0 + 0.5*kRadTolerance/(ySemiAxis*ySemiAxis)) );
290  rad2oi = sqr(p.x()/( xSemiAxis - halfRadTol ))
291  + sqr(p.y()/( ySemiAxis - halfRadTol ));
292 
293  if (rad2oi < sqr( zheight-p.z() ) )
294  {
295  in = ( ( p.z() < -zTopCut + halfRadTol )
296  || ( p.z() > zTopCut - halfRadTol ) ) ? kSurface : kInside;
297  }
298  else
299  {
300  in = kSurface;
301  }
302 
303  return in;
304 }
305 
307 //
308 // Return unit normal of surface closest to p not protected against p=0
309 //
311 {
312 
313  G4double rx = sqr(p.x()/xSemiAxis),
314  ry = sqr(p.y()/ySemiAxis);
315 
316  G4double rds = std::sqrt(rx + ry);
317 
318  G4ThreeVector norm;
319 
320  if( (p.z() < -zTopCut) && ((rx+ry) < sqr(zTopCut + zheight)) )
321  {
322  return G4ThreeVector( 0., 0., -1. );
323  }
324 
325  if( (p.z() > (zheight > zTopCut ? zheight : zTopCut)) &&
326  ((rx+ry) < sqr(zheight-zTopCut)) )
327  {
328  return G4ThreeVector( 0., 0., 1. );
329  }
330 
331  if( p.z() > rds + 2.*zTopCut - zheight )
332  {
333  if ( p.z() > zTopCut )
334  {
335  if( p.x() == 0. )
336  {
337  norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., 1. );
338  return norm /= norm.mag();
339  }
340  if( p.y() == 0. )
341  {
342  norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., 1. );
343  return norm /= norm.mag();
344  }
345 
346  G4double k = std::fabs(p.x()/p.y());
347  G4double c2 = sqr(zheight-zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
348  G4double x = std::sqrt(c2);
349  G4double y = k*x;
350 
351  x /= sqr(xSemiAxis);
352  y /= sqr(ySemiAxis);
353 
354  norm = G4ThreeVector( p.x() < 0. ? -x : x,
355  p.y() < 0. ? -y : y,
356  - ( zheight - zTopCut ) );
357  norm /= norm.mag();
358  norm += G4ThreeVector( 0., 0., 1. );
359  return norm /= norm.mag();
360  }
361 
362  return G4ThreeVector( 0., 0., 1. );
363  }
364 
365  if( p.z() < rds - 2.*zTopCut - zheight )
366  {
367  if( p.x() == 0. )
368  {
369  norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., -1. );
370  return norm /= norm.mag();
371  }
372  if( p.y() == 0. )
373  {
374  norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., -1. );
375  return norm /= norm.mag();
376  }
377 
378  G4double k = std::fabs(p.x()/p.y());
379  G4double c2 = sqr(zheight+zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
380  G4double x = std::sqrt(c2);
381  G4double y = k*x;
382 
383  x /= sqr(xSemiAxis);
384  y /= sqr(ySemiAxis);
385 
386  norm = G4ThreeVector( p.x() < 0. ? -x : x,
387  p.y() < 0. ? -y : y,
388  - ( zheight - zTopCut ) );
389  norm /= norm.mag();
390  norm += G4ThreeVector( 0., 0., -1. );
391  return norm /= norm.mag();
392  }
393 
394  norm = G4ThreeVector(p.x()/sqr(xSemiAxis), p.y()/sqr(ySemiAxis), rds);
395 
396  G4double k = std::tan(pi/8.);
397  G4double c = -zTopCut - k*(zTopCut + zheight);
398 
399  if( p.z() < -k*rds + c )
400  return G4ThreeVector (0.,0.,-1.);
401 
402  return norm /= norm.mag();
403 }
404 
406 //
407 // Calculate distance to shape from outside, along normalised vector
408 // return kInfinity if no intersection, or intersection distance <= tolerance
409 //
411  const G4ThreeVector& v ) const
412 {
413  G4double distMin = kInfinity;
414 
415  // code from EllipticalTube
416 
417  G4double sigz = p.z()+zTopCut;
418 
419  //
420  // Check z = -dz planer surface
421  //
422 
423  if (sigz < halfCarTol)
424  {
425  //
426  // We are "behind" the shape in z, and so can
427  // potentially hit the rear face. Correct direction?
428  //
429  if (v.z() <= 0)
430  {
431  //
432  // As long as we are far enough away, we know we
433  // can't intersect
434  //
435  if (sigz < 0) return kInfinity;
436 
437  //
438  // Otherwise, we don't intersect unless we are
439  // on the surface of the ellipse
440  //
441 
442  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
443  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight+zTopCut ) )
444  return kInfinity;
445 
446  }
447  else
448  {
449  //
450  // How far?
451  //
452  G4double q = -sigz/v.z();
453 
454  //
455  // Where does that place us?
456  //
457  G4double xi = p.x() + q*v.x(),
458  yi = p.y() + q*v.y();
459 
460  //
461  // Is this on the surface (within ellipse)?
462  //
463  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) )
464  {
465  //
466  // Yup. Return q, unless we are on the surface
467  //
468  return (sigz < -halfCarTol) ? q : 0;
469  }
470  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
471  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
472  {
473  //
474  // Else, if we are traveling outwards, we know
475  // we must miss
476  //
477  // return kInfinity;
478  }
479  }
480  }
481 
482  //
483  // Check z = +dz planer surface
484  //
485  sigz = p.z() - zTopCut;
486 
487  if (sigz > -halfCarTol)
488  {
489  if (v.z() >= 0)
490  {
491 
492  if (sigz > 0) return kInfinity;
493 
494  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
495  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight-zTopCut ) )
496  return kInfinity;
497 
498  }
499  else {
500  G4double q = -sigz/v.z();
501 
502  G4double xi = p.x() + q*v.x(),
503  yi = p.y() + q*v.y();
504 
505  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) )
506  {
507  return (sigz > -halfCarTol) ? q : 0;
508  }
509  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
510  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
511  {
512  // return kInfinity;
513  }
514  }
515  }
516 
517 
518 #if 0
519 
520  // check to see if Z plane is relevant
521  //
522  if (p.z() < -zTopCut - 0.5*kCarTolerance)
523  {
524  if (v.z() <= 0.0)
525  return distMin;
526 
527  G4double lambda = (-zTopCut - p.z())/v.z();
528 
529  if ( sqr((lambda*v.x()+p.x())/xSemiAxis) +
530  sqr((lambda*v.y()+p.y())/ySemiAxis) <=
531  sqr(zTopCut + zheight + 0.5*kRadTolerance) )
532  {
533  return distMin = std::fabs(lambda);
534  }
535  }
536 
537  if (p.z() > zTopCut+0.5*kCarTolerance)
538  {
539  if (v.z() >= 0.0)
540  { return distMin; }
541 
542  G4double lambda = (zTopCut - p.z()) / v.z();
543 
544  if ( sqr((lambda*v.x() + p.x())/xSemiAxis) +
545  sqr((lambda*v.y() + p.y())/ySemiAxis) <=
546  sqr(zheight - zTopCut + 0.5*kRadTolerance) )
547  {
548  return distMin = std::fabs(lambda);
549  }
550  }
551 
552  if (p.z() > zTopCut - halfCarTol
553  && p.z() < zTopCut + halfCarTol )
554  {
555  if (v.z() > 0.)
556  { return kInfinity; }
557 
558  return distMin = 0.;
559  }
560 
561  if (p.z() < -zTopCut + halfCarTol
562  && p.z() > -zTopCut - halfCarTol)
563  {
564  if (v.z() < 0.)
565  { return distMin = kInfinity; }
566 
567  return distMin = 0.;
568  }
569 
570 #endif
571 
572  // if we are here then it either intersects or grazes the curved surface
573  // or it does not intersect at all
574  //
575  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
576  G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) +
577  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
578  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) -
579  sqr(zheight - p.z());
580 
581  G4double discr = B*B - 4.*A*C;
582 
583  // if the discriminant is negative it never hits the curved object
584  //
585  if ( discr < -halfCarTol )
586  { return distMin; }
587 
588  // case below is when it hits or grazes the surface
589  //
590  if ( (discr >= - halfCarTol ) && (discr < halfCarTol ) )
591  {
592  return distMin = std::fabs(-B/(2.*A));
593  }
594 
595  G4double plus = (-B+std::sqrt(discr))/(2.*A);
596  G4double minus = (-B-std::sqrt(discr))/(2.*A);
597 
598  // Special case::Point on Surface, Check norm.dot(v)
599 
600  if ( ( std::fabs(plus) < halfCarTol )||( std::fabs(minus) < halfCarTol ) )
601  {
602  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
603  p.y()/(ySemiAxis*ySemiAxis),
604  -( p.z() - zheight ));
605  if ( truenorm*v >= 0) // going outside the solid from surface
606  {
607  return kInfinity;
608  }
609  else
610  {
611  return 0;
612  }
613  }
614 
615  // G4double lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
616  G4double lambda = 0;
617 
618  if ( minus > halfCarTol && minus < distMin )
619  {
620  lambda = minus ;
621  // check normal vector n * v < 0
622  G4ThreeVector pin = p + lambda*v;
623  if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
624  {
625  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
626  pin.y()/(ySemiAxis*ySemiAxis),
627  - ( pin.z() - zheight ));
628  if ( truenorm*v < 0)
629  { // yes, going inside the solid
630  distMin = lambda;
631  }
632  }
633  }
634  if ( plus > halfCarTol && plus < distMin )
635  {
636  lambda = plus ;
637  // check normal vector n * v < 0
638  G4ThreeVector pin = p + lambda*v;
639  if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
640  {
641  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
642  pin.y()/(ySemiAxis*ySemiAxis),
643  - ( pin.z() - zheight ) );
644  if ( truenorm*v < 0)
645  { // yes, going inside the solid
646  distMin = lambda;
647  }
648  }
649  }
650  if (distMin < halfCarTol) distMin=0.;
651  return distMin ;
652 }
653 
655 //
656 // Calculate distance (<= actual) to closest surface of shape from outside
657 // Return 0 if point inside
658 //
660 {
661  G4double distR, distR2, distZ, maxDim;
662  G4double distRad;
663 
664  // check if the point lies either below z=-zTopCut in bottom elliptical
665  // region or on top within cut elliptical region
666  //
667  if( (p.z() <= -zTopCut) && (sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
668  <= sqr(zTopCut + zheight + 0.5*kCarTolerance )) )
669  {
670  //return distZ = std::fabs(zTopCut - p.z());
671  return distZ = std::fabs(zTopCut + p.z());
672  }
673 
674  if( (p.z() >= zTopCut) && (sqr(p.x()/xSemiAxis)+sqr(p.y()/ySemiAxis)
675  <= sqr(zheight - zTopCut + kCarTolerance/2.0 )) )
676  {
677  return distZ = std::fabs(p.z() - zTopCut);
678  }
679 
680  // below we use the following approximation: we take the largest of the
681  // axes and find the shortest distance to the circular (cut) cone of that
682  // radius.
683  //
684  maxDim = xSemiAxis >= ySemiAxis ? xSemiAxis:ySemiAxis;
685  distRad = std::sqrt(p.x()*p.x()+p.y()*p.y());
686 
687  if( p.z() > maxDim*distRad + zTopCut*(1.+maxDim)-sqr(maxDim)*zheight )
688  {
689  distR2 = sqr(p.z() - zTopCut) + sqr(distRad - maxDim*(zheight - zTopCut));
690  return std::sqrt( distR2 );
691  }
692 
693  if( distRad > maxDim*( zheight - p.z() ) )
694  {
695  if( p.z() > maxDim*distRad - (zTopCut*(1.+maxDim)+sqr(maxDim)*zheight) )
696  {
697  G4double zVal = (p.z()-maxDim*(distRad-maxDim*zheight))/(1.+sqr(maxDim));
698  G4double rVal = maxDim*(zheight - zVal);
699  return distR = std::sqrt(sqr(p.z() - zVal) + sqr(distRad - rVal));
700  }
701  }
702 
703  if( distRad <= maxDim*(zheight - p.z()) )
704  {
705  distR2 = sqr(distRad - maxDim*(zheight + zTopCut)) + sqr(p.z() + zTopCut);
706  return std::sqrt( distR2 );
707  }
708 
709  return distR = 0;
710 }
711 
713 //
714 // Calculate distance to surface of shape from `inside',
715 // allowing for tolerance
716 //
718  const G4ThreeVector& v,
719  const G4bool calcNorm,
720  G4bool *validNorm,
721  G4ThreeVector *n ) const
722 {
723  G4double distMin, lambda;
724  enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
725 
726  distMin = kInfinity;
727  surface = kNoSurf;
728 
729  if (v.z() < 0.0)
730  {
731  lambda = (-p.z() - zTopCut)/v.z();
732 
733  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) +
734  sqr((p.y() + lambda*v.y())/ySemiAxis)) <
735  sqr(zheight + zTopCut + 0.5*kCarTolerance) )
736  {
737  distMin = std::fabs(lambda);
738 
739  if (!calcNorm) { return distMin; }
740  }
741  distMin = std::fabs(lambda);
742  surface = kPlaneSurf;
743  }
744 
745  if (v.z() > 0.0)
746  {
747  lambda = (zTopCut - p.z()) / v.z();
748 
749  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis)
750  + sqr((p.y() + lambda*v.y())/ySemiAxis) )
751  < (sqr(zheight - zTopCut + 0.5*kCarTolerance)) )
752  {
753  distMin = std::fabs(lambda);
754  if (!calcNorm) { return distMin; }
755  }
756  distMin = std::fabs(lambda);
757  surface = kPlaneSurf;
758  }
759 
760  // if we are here then it either intersects or grazes the
761  // curved surface...
762  //
763  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
764  G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) +
765  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
766  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
767  - sqr(zheight - p.z());
768 
769  G4double discr = B*B - 4.*A*C;
770 
771  if ( discr >= - 0.5*kCarTolerance && discr < 0.5*kCarTolerance )
772  {
773  if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); }
774  }
775 
776  else if ( discr > 0.5*kCarTolerance )
777  {
778  G4double plus = (-B+std::sqrt(discr))/(2.*A);
779  G4double minus = (-B-std::sqrt(discr))/(2.*A);
780 
781  if ( plus > 0.5*kCarTolerance && minus > 0.5*kCarTolerance )
782  {
783  // take the shorter distance
784  //
785  lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
786  }
787  else
788  {
789  // at least one solution is close to zero or negative
790  // so, take small positive solution or zero
791  //
792  lambda = plus > -0.5*kCarTolerance ? plus : 0;
793  }
794 
795  if ( std::fabs(lambda) < distMin )
796  {
797  if( std::fabs(lambda) > 0.5*kCarTolerance)
798  {
799  distMin = std::fabs(lambda);
800  surface = kCurvedSurf;
801  }
802  else // Point is On the Surface, Check Normal
803  {
804  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
805  p.y()/(ySemiAxis*ySemiAxis),
806  -( p.z() - zheight ));
807  if( truenorm.dot(v) > 0 )
808  {
809  distMin = 0.0;
810  surface = kCurvedSurf;
811  }
812  }
813  }
814  }
815 
816  // set normal if requested
817  //
818  if (calcNorm)
819  {
820  if (surface == kNoSurf)
821  {
822  *validNorm = false;
823  }
824  else
825  {
826  *validNorm = true;
827  switch (surface)
828  {
829  case kPlaneSurf:
830  {
831  *n = G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
832  }
833  break;
834 
835  case kCurvedSurf:
836  {
837  G4ThreeVector pexit = p + distMin*v;
838  G4ThreeVector truenorm( pexit.x()/(xSemiAxis*xSemiAxis),
839  pexit.y()/(ySemiAxis*ySemiAxis),
840  -( pexit.z() - zheight ) );
841  truenorm /= truenorm.mag();
842  *n= truenorm;
843  }
844  break;
845 
846  default: // Should never reach this case ...
847  DumpInfo();
848  std::ostringstream message;
849  G4int oldprc = message.precision(16);
850  message << "Undefined side for valid surface normal to solid."
851  << G4endl
852  << "Position:" << G4endl
853  << " p.x() = " << p.x()/mm << " mm" << G4endl
854  << " p.y() = " << p.y()/mm << " mm" << G4endl
855  << " p.z() = " << p.z()/mm << " mm" << G4endl
856  << "Direction:" << G4endl
857  << " v.x() = " << v.x() << G4endl
858  << " v.y() = " << v.y() << G4endl
859  << " v.z() = " << v.z() << G4endl
860  << "Proposed distance :" << G4endl
861  << " distMin = " << distMin/mm << " mm";
862  message.precision(oldprc);
863  G4Exception("G4EllipticalCone::DistanceToOut(p,v,..)",
864  "GeomSolids1002", JustWarning, message);
865  break;
866  }
867  }
868  }
869 
870  if (distMin<0.5*kCarTolerance) { distMin=0; }
871 
872  return distMin;
873 }
874 
876 //
877 // Calculate distance (<=actual) to closest surface of shape from inside
878 //
880 {
881 #ifdef G4SPECSDEBUG
882  if( Inside(p) == kOutside )
883  {
884  DumpInfo();
885  std::ostringstream message;
886  G4int oldprc = message.precision(16);
887  message << "Point p is outside !?" << G4endl
888  << "Position:" << G4endl
889  << " p.x() = " << p.x()/mm << " mm" << G4endl
890  << " p.y() = " << p.y()/mm << " mm" << G4endl
891  << " p.z() = " << p.z()/mm << " mm";
892  message.precision(oldprc) ;
893  G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
894  JustWarning, message);
895  }
896 #endif
897  // The safety is calculated in the scaled space where elliptical cone
898  // becomes a circular cone with radius equal to the smaller of the axes
899  //
900  G4double px = p.x(), py = p.y(), pz = p.z();
901  G4double axis;
902  if (xSemiAxis < ySemiAxis)
903  {
904  axis = xSemiAxis;
905  py *= xSemiAxis/ySemiAxis; // scale y
906  }
907  else
908  {
909  axis = ySemiAxis;
910  px *= ySemiAxis/xSemiAxis; // scale x
911  }
912 
913  G4double distZ = zTopCut - std::abs(pz) ;
914  if (distZ <= 0) return 0; // point is outside
915 
916  G4double rho = axis*(zheight-pz); // radius at z = p.z()
917  G4double pr = std::sqrt(px*px+py*py);
918  if (pr >= rho) return 0; // point is outside
919 
920  G4double distR = (rho-pr)/std::sqrt(1+axis*axis);
921  return std::min(distR,distZ);
922 }
923 
925 //
926 // GetEntityType
927 //
929 {
930  return G4String("G4EllipticalCone");
931 }
932 
934 //
935 // Make a clone of the object
936 //
938 {
939  return new G4EllipticalCone(*this);
940 }
941 
943 //
944 // Stream object contents to an output stream
945 //
946 std::ostream& G4EllipticalCone::StreamInfo( std::ostream& os ) const
947 {
948  G4int oldprc = os.precision(16);
949  os << "-----------------------------------------------------------\n"
950  << " *** Dump for solid - " << GetName() << " ***\n"
951  << " ===================================================\n"
952  << " Solid type: G4EllipticalCone\n"
953  << " Parameters: \n"
954 
955  << " semi-axis x: " << xSemiAxis/mm << " mm \n"
956  << " semi-axis y: " << ySemiAxis/mm << " mm \n"
957  << " height z: " << zheight/mm << " mm \n"
958  << " half length in z: " << zTopCut/mm << " mm \n"
959  << "-----------------------------------------------------------\n";
960  os.precision(oldprc);
961 
962  return os;
963 }
964 
966 //
967 // GetPointOnSurface
968 //
969 // returns quasi-uniformly distributed point on surface of elliptical cone
970 //
972 {
973 
974  G4double phi, sinphi, cosphi, aOne, aTwo, aThree,
975  chose, zRand, rRand1, rRand2;
976 
977  G4double rOne = std::sqrt(sqr(xSemiAxis)
978  + sqr(ySemiAxis))*(zheight - zTopCut);
979  G4double rTwo = std::sqrt(sqr(xSemiAxis)
980  + sqr(ySemiAxis))*(zheight + zTopCut);
981 
982  G4int it1=0, it2=0;
983 
984  aOne = pi*(rOne + rTwo)*std::sqrt(sqr(rOne - rTwo)+sqr(2.*zTopCut));
985  aTwo = pi*xSemiAxis*ySemiAxis*sqr(zheight+zTopCut);
986  aThree = pi*xSemiAxis*ySemiAxis*sqr(zheight-zTopCut);
987 
988  phi = G4RandFlat::shoot(0.,twopi);
989  cosphi = std::cos(phi);
990  sinphi = std::sin(phi);
991 
992  if(zTopCut >= zheight) aThree = 0.;
993 
994  chose = G4RandFlat::shoot(0.,aOne+aTwo+aThree);
995  if((chose>=0.) && (chose<aOne))
996  {
997  zRand = G4RandFlat::shoot(-zTopCut,zTopCut);
998  return G4ThreeVector(xSemiAxis*(zheight-zRand)*cosphi,
999  ySemiAxis*(zheight-zRand)*sinphi,zRand);
1000  }
1001  else if((chose>=aOne) && (chose<aOne+aTwo))
1002  {
1003  do // Loop checking, 13.08.2015, G.Cosmo
1004  {
1005  rRand1 = G4RandFlat::shoot(0.,1.) ;
1006  rRand2 = G4RandFlat::shoot(0.,1.) ;
1007  } while (( rRand2 >= rRand1 ) && (++it1 < 1000)) ;
1008 
1009  return G4ThreeVector(rRand1*xSemiAxis*(zheight+zTopCut)*cosphi,
1010  rRand1*ySemiAxis*(zheight+zTopCut)*sinphi, -zTopCut);
1011 
1012  }
1013  // else
1014  //
1015 
1016  do // Loop checking, 13.08.2015, G.Cosmo
1017  {
1018  rRand1 = G4RandFlat::shoot(0.,1.) ;
1019  rRand2 = G4RandFlat::shoot(0.,1.) ;
1020  } while (( rRand2 >= rRand1 ) && (++it2 < 1000));
1021 
1022  return G4ThreeVector(rRand1*xSemiAxis*(zheight-zTopCut)*cosphi,
1023  rRand1*ySemiAxis*(zheight-zTopCut)*sinphi, zTopCut);
1024 }
1025 
1026 //
1027 // Methods for visualisation
1028 //
1029 
1031 {
1032  scene.AddSolid(*this);
1033 }
1034 
1036 {
1037  // Define the sides of the box into which the solid instance would fit.
1038  //
1039  G4ThreeVector pmin,pmax;
1040  Extent(pmin,pmax);
1041  return G4VisExtent(pmin.x(),pmax.x(),
1042  pmin.y(),pmax.y(),
1043  pmin.z(),pmax.z());
1044 }
1045 
1047 {
1048  return new G4PolyhedronEllipticalCone(xSemiAxis, ySemiAxis, zheight, zTopCut);
1049 }
1050 
1052 {
1053  if ( (!fpPolyhedron)
1057  {
1058  G4AutoLock l(&polyhedronMutex);
1059  delete fpPolyhedron;
1061  fRebuildPolyhedron = false;
1062  l.unlock();
1063  }
1064  return fpPolyhedron;
1065 }
void set(double x, double y, double z)
G4String GetName() const
G4VSolid * Clone() const
ThreeVector shoot(const G4int Ap, const G4int Af)
static constexpr double mm
Definition: G4SIunits.hh:115
G4EllipticalCone & operator=(const G4EllipticalCone &rhs)
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetSemiAxisX() const
CLHEP::Hep3Vector G4ThreeVector
double x() const
double dot(const Hep3Vector &) const
const char * p
Definition: xmltok.h:285
EInside Inside(const G4ThreeVector &p) const
double C(double temp)
double B(double temperature)
virtual void AddSolid(const G4Box &)=0
G4Polyhedron * CreatePolyhedron() const
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
double z() const
void DumpInfo() const
static constexpr double twopi
Definition: G4SIunits.hh:76
std::ostream & StreamInfo(std::ostream &os) const
G4ThreeVector GetPointOnSurface() const
double A(double temperature)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4Polyhedron * fpPolyhedron
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4EllipticalCone(const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double zMax, G4double pzTopCut)
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4double GetZMax() const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
std::vector< G4ThreeVector > G4ThreeVectorList
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4double GetSemiAxisY() const
G4int G4Mutex
Definition: G4Threading.hh:173
void SetSemiAxis(G4double x, G4double y, G4double z)
static G4int GetNumberOfRotationSteps()
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
void SetZCut(G4double newzTopCut)
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetZTopCut() const
#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
double G4double
Definition: G4Types.hh:76
G4Polyhedron * GetPolyhedron() const
G4GeometryType GetEntityType() const
G4VisExtent GetExtent() const
double mag() const
virtual ~G4EllipticalCone()
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
static G4GeometryTolerance * GetInstance()