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