Geant4  10.02.p01
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 92024 2015-08-13 14:16:00Z 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 "G4SolidExtentList.hh"
46 #include "G4VoxelLimits.hh"
47 #include "G4AffineTransform.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 
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  //
161  xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
163  fRebuildPolyhedron = false;
164  delete fpPolyhedron; fpPolyhedron = 0;
165 
166  return *this;
167 }
168 
170 //
171 // Calculate extent under transform and specified limit
172 //
173 G4bool
175  const G4VoxelLimits &voxelLimit,
176  const G4AffineTransform &transform,
177  G4double &min, G4double &max ) const
178 {
179  G4SolidExtentList extentList( axis, voxelLimit );
180 
181  //
182  // We are going to divide up our elliptical face into small pieces
183  //
184 
185  //
186  // Choose phi size of our segment(s) based on constants as
187  // defined in meshdefs.hh
188  //
189  G4int numPhi = kMaxMeshSections;
190  G4double sigPhi = twopi/numPhi;
191 
192  //
193  // We have to be careful to keep our segments completely outside
194  // of the elliptical surface. To do so we imagine we have
195  // a simple (unit radius) circular cross section (as in G4Tubs)
196  // and then "stretch" the dimensions as necessary to fit the ellipse.
197  //
198  G4double rFudge = 1.0/std::cos(0.5*sigPhi);
199  G4double dxFudgeBot = xSemiAxis*2.*zheight*rFudge,
200  dyFudgeBot = ySemiAxis*2.*zheight*rFudge;
201  G4double dxFudgeTop = xSemiAxis*(zheight-zTopCut)*rFudge,
202  dyFudgeTop = ySemiAxis*(zheight-zTopCut)*rFudge;
203 
204  //
205  // As we work around the elliptical surface, we build
206  // a "phi" segment on the way, and keep track of two
207  // additional polygons for the two ends.
208  //
209  G4ClippablePolygon endPoly1, endPoly2, phiPoly;
210 
211  G4double phi = 0,
212  cosPhi = std::cos(phi),
213  sinPhi = std::sin(phi);
214  G4ThreeVector v0( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut ),
215  v1( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut ),
216  w0, w1;
217  transform.ApplyPointTransform( v0 );
218  transform.ApplyPointTransform( v1 );
219  do // Loop checking, 13.08.2015, G.Cosmo
220  {
221  phi += sigPhi;
222  if (numPhi == 1) phi = 0; // Try to avoid roundoff
223  cosPhi = std::cos(phi),
224  sinPhi = std::sin(phi);
225 
226  w0 = G4ThreeVector( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut );
227  w1 = G4ThreeVector( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut );
228  transform.ApplyPointTransform( w0 );
229  transform.ApplyPointTransform( w1 );
230 
231  //
232  // Add a point to our z ends
233  //
234  endPoly1.AddVertexInOrder( v0 );
235  endPoly2.AddVertexInOrder( v1 );
236 
237  //
238  // Build phi polygon
239  //
240  phiPoly.ClearAllVertices();
241 
242  phiPoly.AddVertexInOrder( v0 );
243  phiPoly.AddVertexInOrder( v1 );
244  phiPoly.AddVertexInOrder( w1 );
245  phiPoly.AddVertexInOrder( w0 );
246 
247  if (phiPoly.PartialClip( voxelLimit, axis ))
248  {
249  //
250  // Get unit normal
251  //
252  phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
253 
254  extentList.AddSurface( phiPoly );
255  }
256 
257  //
258  // Next vertex
259  //
260  v0 = w0;
261  v1 = w1;
262  } while( --numPhi > 0 );
263 
264  //
265  // Process the end pieces
266  //
267  if (endPoly1.PartialClip( voxelLimit, axis ))
268  {
269  static const G4ThreeVector normal(0,0,+1);
270  endPoly1.SetNormal( transform.TransformAxis(normal) );
271  extentList.AddSurface( endPoly1 );
272  }
273 
274  if (endPoly2.PartialClip( voxelLimit, axis ))
275  {
276  static const G4ThreeVector normal(0,0,-1);
277  endPoly2.SetNormal( transform.TransformAxis(normal) );
278  extentList.AddSurface( endPoly2 );
279  }
280 
281  //
282  // Return min/max value
283  //
284  return extentList.GetExtent( min, max );
285 }
286 
288 //
289 // Return whether point inside/outside/on surface
290 // Split into radius, phi, theta checks
291 // Each check modifies `in', or returns as approprate
292 //
294 {
295  G4double rad2oo, // outside surface outer tolerance
296  rad2oi; // outside surface inner tolerance
297 
298  EInside in;
299 
300  // check this side of z cut first, because that's fast
301  //
302 
303  if ( (p.z() < -zTopCut - halfCarTol)
304  || (p.z() > zTopCut + halfCarTol ) )
305  {
306  return in = kOutside;
307  }
308 
309  rad2oo= sqr(p.x()/( xSemiAxis + halfRadTol ))
310  + sqr(p.y()/( ySemiAxis + halfRadTol ));
311 
312  if ( rad2oo > sqr( zheight-p.z() ) )
313  {
314  return in = kOutside;
315  }
316 
317  // rad2oi= sqr( p.x()*(1.0 + 0.5*kRadTolerance/(xSemiAxis*xSemiAxis)) )
318  // + sqr( p.y()*(1.0 + 0.5*kRadTolerance/(ySemiAxis*ySemiAxis)) );
319  rad2oi = sqr(p.x()/( xSemiAxis - halfRadTol ))
320  + sqr(p.y()/( ySemiAxis - halfRadTol ));
321 
322  if (rad2oi < sqr( zheight-p.z() ) )
323  {
324  in = ( ( p.z() < -zTopCut + halfRadTol )
325  || ( p.z() > zTopCut - halfRadTol ) ) ? kSurface : kInside;
326  }
327  else
328  {
329  in = kSurface;
330  }
331 
332  return in;
333 }
334 
336 //
337 // Return unit normal of surface closest to p not protected against p=0
338 //
340 {
341 
342  G4double rx = sqr(p.x()/xSemiAxis),
343  ry = sqr(p.y()/ySemiAxis);
344 
345  G4double rds = std::sqrt(rx + ry);
346 
347  G4ThreeVector norm;
348 
349  if( (p.z() < -zTopCut) && ((rx+ry) < sqr(zTopCut + zheight)) )
350  {
351  return G4ThreeVector( 0., 0., -1. );
352  }
353 
354  if( (p.z() > (zheight > zTopCut ? zheight : zTopCut)) &&
355  ((rx+ry) < sqr(zheight-zTopCut)) )
356  {
357  return G4ThreeVector( 0., 0., 1. );
358  }
359 
360  if( p.z() > rds + 2.*zTopCut - zheight )
361  {
362  if ( p.z() > zTopCut )
363  {
364  if( p.x() == 0. )
365  {
366  norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., 1. );
367  return norm /= norm.mag();
368  }
369  if( p.y() == 0. )
370  {
371  norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., 1. );
372  return norm /= norm.mag();
373  }
374 
375  G4double k = std::fabs(p.x()/p.y());
377  G4double x = std::sqrt(c2);
378  G4double y = k*x;
379 
380  x /= sqr(xSemiAxis);
381  y /= sqr(ySemiAxis);
382 
383  norm = G4ThreeVector( p.x() < 0. ? -x : x,
384  p.y() < 0. ? -y : y,
385  - ( zheight - zTopCut ) );
386  norm /= norm.mag();
387  norm += G4ThreeVector( 0., 0., 1. );
388  return norm /= norm.mag();
389  }
390 
391  return G4ThreeVector( 0., 0., 1. );
392  }
393 
394  if( p.z() < rds - 2.*zTopCut - zheight )
395  {
396  if( p.x() == 0. )
397  {
398  norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., -1. );
399  return norm /= norm.mag();
400  }
401  if( p.y() == 0. )
402  {
403  norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., -1. );
404  return norm /= norm.mag();
405  }
406 
407  G4double k = std::fabs(p.x()/p.y());
409  G4double x = std::sqrt(c2);
410  G4double y = k*x;
411 
412  x /= sqr(xSemiAxis);
413  y /= sqr(ySemiAxis);
414 
415  norm = G4ThreeVector( p.x() < 0. ? -x : x,
416  p.y() < 0. ? -y : y,
417  - ( zheight - zTopCut ) );
418  norm /= norm.mag();
419  norm += G4ThreeVector( 0., 0., -1. );
420  return norm /= norm.mag();
421  }
422 
423  norm = G4ThreeVector(p.x()/sqr(xSemiAxis), p.y()/sqr(ySemiAxis), rds);
424 
425  G4double k = std::tan(pi/8.);
426  G4double c = -zTopCut - k*(zTopCut + zheight);
427 
428  if( p.z() < -k*rds + c )
429  return G4ThreeVector (0.,0.,-1.);
430 
431  return norm /= norm.mag();
432 }
433 
435 //
436 // Calculate distance to shape from outside, along normalised vector
437 // return kInfinity if no intersection, or intersection distance <= tolerance
438 //
440  const G4ThreeVector& v ) const
441 {
442  G4double distMin = kInfinity;
443 
444  // code from EllipticalTube
445 
446  G4double sigz = p.z()+zTopCut;
447 
448  //
449  // Check z = -dz planer surface
450  //
451 
452  if (sigz < halfCarTol)
453  {
454  //
455  // We are "behind" the shape in z, and so can
456  // potentially hit the rear face. Correct direction?
457  //
458  if (v.z() <= 0)
459  {
460  //
461  // As long as we are far enough away, we know we
462  // can't intersect
463  //
464  if (sigz < 0) return kInfinity;
465 
466  //
467  // Otherwise, we don't intersect unless we are
468  // on the surface of the ellipse
469  //
470 
471  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
472  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight+zTopCut ) )
473  return kInfinity;
474 
475  }
476  else
477  {
478  //
479  // How far?
480  //
481  G4double q = -sigz/v.z();
482 
483  //
484  // Where does that place us?
485  //
486  G4double xi = p.x() + q*v.x(),
487  yi = p.y() + q*v.y();
488 
489  //
490  // Is this on the surface (within ellipse)?
491  //
492  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) )
493  {
494  //
495  // Yup. Return q, unless we are on the surface
496  //
497  return (sigz < -halfCarTol) ? q : 0;
498  }
499  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
500  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
501  {
502  //
503  // Else, if we are traveling outwards, we know
504  // we must miss
505  //
506  // return kInfinity;
507  }
508  }
509  }
510 
511  //
512  // Check z = +dz planer surface
513  //
514  sigz = p.z() - zTopCut;
515 
516  if (sigz > -halfCarTol)
517  {
518  if (v.z() >= 0)
519  {
520 
521  if (sigz > 0) return kInfinity;
522 
523  if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
524  + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight-zTopCut ) )
525  return kInfinity;
526 
527  }
528  else {
529  G4double q = -sigz/v.z();
530 
531  G4double xi = p.x() + q*v.x(),
532  yi = p.y() + q*v.y();
533 
534  if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) )
535  {
536  return (sigz > -halfCarTol) ? q : 0;
537  }
538  else if (xi/(xSemiAxis*xSemiAxis)*v.x()
539  + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
540  {
541  // return kInfinity;
542  }
543  }
544  }
545 
546 
547 #if 0
548 
549  // check to see if Z plane is relevant
550  //
551  if (p.z() < -zTopCut - 0.5*kCarTolerance)
552  {
553  if (v.z() <= 0.0)
554  return distMin;
555 
556  G4double lambda = (-zTopCut - p.z())/v.z();
557 
558  if ( sqr((lambda*v.x()+p.x())/xSemiAxis) +
559  sqr((lambda*v.y()+p.y())/ySemiAxis) <=
560  sqr(zTopCut + zheight + 0.5*kRadTolerance) )
561  {
562  return distMin = std::fabs(lambda);
563  }
564  }
565 
566  if (p.z() > zTopCut+0.5*kCarTolerance)
567  {
568  if (v.z() >= 0.0)
569  { return distMin; }
570 
571  G4double lambda = (zTopCut - p.z()) / v.z();
572 
573  if ( sqr((lambda*v.x() + p.x())/xSemiAxis) +
574  sqr((lambda*v.y() + p.y())/ySemiAxis) <=
575  sqr(zheight - zTopCut + 0.5*kRadTolerance) )
576  {
577  return distMin = std::fabs(lambda);
578  }
579  }
580 
581  if (p.z() > zTopCut - halfCarTol
582  && p.z() < zTopCut + halfCarTol )
583  {
584  if (v.z() > 0.)
585  { return kInfinity; }
586 
587  return distMin = 0.;
588  }
589 
590  if (p.z() < -zTopCut + halfCarTol
591  && p.z() > -zTopCut - halfCarTol)
592  {
593  if (v.z() < 0.)
594  { return distMin = kInfinity; }
595 
596  return distMin = 0.;
597  }
598 
599 #endif
600 
601  // if we are here then it either intersects or grazes the curved surface
602  // or it does not intersect at all
603  //
604  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
605  G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) +
606  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
607  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) -
608  sqr(zheight - p.z());
609 
610  G4double discr = B*B - 4.*A*C;
611 
612  // if the discriminant is negative it never hits the curved object
613  //
614  if ( discr < -halfCarTol )
615  { return distMin; }
616 
617  // case below is when it hits or grazes the surface
618  //
619  if ( (discr >= - halfCarTol ) && (discr < halfCarTol ) )
620  {
621  return distMin = std::fabs(-B/(2.*A));
622  }
623 
624  G4double plus = (-B+std::sqrt(discr))/(2.*A);
625  G4double minus = (-B-std::sqrt(discr))/(2.*A);
626 
627  // Special case::Point on Surface, Check norm.dot(v)
628 
629  if ( ( std::fabs(plus) < halfCarTol )||( std::fabs(minus) < halfCarTol ) )
630  {
631  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
632  p.y()/(ySemiAxis*ySemiAxis),
633  -( p.z() - zheight ));
634  if ( truenorm*v >= 0) // going outside the solid from surface
635  {
636  return kInfinity;
637  }
638  else
639  {
640  return 0;
641  }
642  }
643 
644  // G4double lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
645  G4double lambda = 0;
646 
647  if ( minus > halfCarTol && minus < distMin )
648  {
649  lambda = minus ;
650  // check normal vector n * v < 0
651  G4ThreeVector pin = p + lambda*v;
652  if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
653  {
654  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
655  pin.y()/(ySemiAxis*ySemiAxis),
656  - ( pin.z() - zheight ));
657  if ( truenorm*v < 0)
658  { // yes, going inside the solid
659  distMin = lambda;
660  }
661  }
662  }
663  if ( plus > halfCarTol && plus < distMin )
664  {
665  lambda = plus ;
666  // check normal vector n * v < 0
667  G4ThreeVector pin = p + lambda*v;
668  if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
669  {
670  G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
671  pin.y()/(ySemiAxis*ySemiAxis),
672  - ( pin.z() - zheight ) );
673  if ( truenorm*v < 0)
674  { // yes, going inside the solid
675  distMin = lambda;
676  }
677  }
678  }
679  if (distMin < halfCarTol) distMin=0.;
680  return distMin ;
681 }
682 
684 //
685 // Calculate distance (<= actual) to closest surface of shape from outside
686 // Return 0 if point inside
687 //
689 {
690  G4double distR, distR2, distZ, maxDim;
691  G4double distRad;
692 
693  // check if the point lies either below z=-zTopCut in bottom elliptical
694  // region or on top within cut elliptical region
695  //
696  if( (p.z() <= -zTopCut) && (sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
697  <= sqr(zTopCut + zheight + 0.5*kCarTolerance )) )
698  {
699  //return distZ = std::fabs(zTopCut - p.z());
700  return distZ = std::fabs(zTopCut + p.z());
701  }
702 
703  if( (p.z() >= zTopCut) && (sqr(p.x()/xSemiAxis)+sqr(p.y()/ySemiAxis)
704  <= sqr(zheight - zTopCut + kCarTolerance/2.0 )) )
705  {
706  return distZ = std::fabs(p.z() - zTopCut);
707  }
708 
709  // below we use the following approximation: we take the largest of the
710  // axes and find the shortest distance to the circular (cut) cone of that
711  // radius.
712  //
713  maxDim = xSemiAxis >= ySemiAxis ? xSemiAxis:ySemiAxis;
714  distRad = std::sqrt(p.x()*p.x()+p.y()*p.y());
715 
716  if( p.z() > maxDim*distRad + zTopCut*(1.+maxDim)-sqr(maxDim)*zheight )
717  {
718  distR2 = sqr(p.z() - zTopCut) + sqr(distRad - maxDim*(zheight - zTopCut));
719  return std::sqrt( distR2 );
720  }
721 
722  if( distRad > maxDim*( zheight - p.z() ) )
723  {
724  if( p.z() > maxDim*distRad - (zTopCut*(1.+maxDim)+sqr(maxDim)*zheight) )
725  {
726  G4double zVal = (p.z()-maxDim*(distRad-maxDim*zheight))/(1.+sqr(maxDim));
727  G4double rVal = maxDim*(zheight - zVal);
728  return distR = std::sqrt(sqr(p.z() - zVal) + sqr(distRad - rVal));
729  }
730  }
731 
732  if( distRad <= maxDim*(zheight - p.z()) )
733  {
734  distR2 = sqr(distRad - maxDim*(zheight + zTopCut)) + sqr(p.z() + zTopCut);
735  return std::sqrt( distR2 );
736  }
737 
738  return distR = 0;
739 }
740 
742 //
743 // Calculate distance to surface of shape from `inside',
744 // allowing for tolerance
745 //
747  const G4ThreeVector& v,
748  const G4bool calcNorm,
749  G4bool *validNorm,
750  G4ThreeVector *n ) const
751 {
752  G4double distMin, lambda;
753  enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
754 
755  distMin = kInfinity;
756  surface = kNoSurf;
757 
758  if (v.z() < 0.0)
759  {
760  lambda = (-p.z() - zTopCut)/v.z();
761 
762  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) +
763  sqr((p.y() + lambda*v.y())/ySemiAxis)) <
764  sqr(zheight + zTopCut + 0.5*kCarTolerance) )
765  {
766  distMin = std::fabs(lambda);
767 
768  if (!calcNorm) { return distMin; }
769  }
770  distMin = std::fabs(lambda);
771  surface = kPlaneSurf;
772  }
773 
774  if (v.z() > 0.0)
775  {
776  lambda = (zTopCut - p.z()) / v.z();
777 
778  if ( (sqr((p.x() + lambda*v.x())/xSemiAxis)
779  + sqr((p.y() + lambda*v.y())/ySemiAxis) )
780  < (sqr(zheight - zTopCut + 0.5*kCarTolerance)) )
781  {
782  distMin = std::fabs(lambda);
783  if (!calcNorm) { return distMin; }
784  }
785  distMin = std::fabs(lambda);
786  surface = kPlaneSurf;
787  }
788 
789  // if we are here then it either intersects or grazes the
790  // curved surface...
791  //
792  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
793  G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) +
794  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
795  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
796  - sqr(zheight - p.z());
797 
798  G4double discr = B*B - 4.*A*C;
799 
800  if ( discr >= - 0.5*kCarTolerance && discr < 0.5*kCarTolerance )
801  {
802  if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); }
803  }
804 
805  else if ( discr > 0.5*kCarTolerance )
806  {
807  G4double plus = (-B+std::sqrt(discr))/(2.*A);
808  G4double minus = (-B-std::sqrt(discr))/(2.*A);
809 
810  if ( plus > 0.5*kCarTolerance && minus > 0.5*kCarTolerance )
811  {
812  // take the shorter distance
813  //
814  lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
815  }
816  else
817  {
818  // at least one solution is close to zero or negative
819  // so, take small positive solution or zero
820  //
821  lambda = plus > -0.5*kCarTolerance ? plus : 0;
822  }
823 
824  if ( std::fabs(lambda) < distMin )
825  {
826  if( std::fabs(lambda) > 0.5*kCarTolerance)
827  {
828  distMin = std::fabs(lambda);
829  surface = kCurvedSurf;
830  }
831  else // Point is On the Surface, Check Normal
832  {
833  G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
834  p.y()/(ySemiAxis*ySemiAxis),
835  -( p.z() - zheight ));
836  if( truenorm.dot(v) > 0 )
837  {
838  distMin = 0.0;
839  surface = kCurvedSurf;
840  }
841  }
842  }
843  }
844 
845  // set normal if requested
846  //
847  if (calcNorm)
848  {
849  if (surface == kNoSurf)
850  {
851  *validNorm = false;
852  }
853  else
854  {
855  *validNorm = true;
856  switch (surface)
857  {
858  case kPlaneSurf:
859  {
860  *n = G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
861  }
862  break;
863 
864  case kCurvedSurf:
865  {
866  G4ThreeVector pexit = p + distMin*v;
867  G4ThreeVector truenorm( pexit.x()/(xSemiAxis*xSemiAxis),
868  pexit.y()/(ySemiAxis*ySemiAxis),
869  -( pexit.z() - zheight ) );
870  truenorm /= truenorm.mag();
871  *n= truenorm;
872  }
873  break;
874 
875  default: // Should never reach this case ...
876  DumpInfo();
877  std::ostringstream message;
878  G4int oldprc = message.precision(16);
879  message << "Undefined side for valid surface normal to solid."
880  << G4endl
881  << "Position:" << G4endl
882  << " p.x() = " << p.x()/mm << " mm" << G4endl
883  << " p.y() = " << p.y()/mm << " mm" << G4endl
884  << " p.z() = " << p.z()/mm << " mm" << G4endl
885  << "Direction:" << G4endl
886  << " v.x() = " << v.x() << G4endl
887  << " v.y() = " << v.y() << G4endl
888  << " v.z() = " << v.z() << G4endl
889  << "Proposed distance :" << G4endl
890  << " distMin = " << distMin/mm << " mm";
891  message.precision(oldprc);
892  G4Exception("G4EllipticalCone::DistanceToOut(p,v,..)",
893  "GeomSolids1002", JustWarning, message);
894  break;
895  }
896  }
897  }
898 
899  if (distMin<0.5*kCarTolerance) { distMin=0; }
900 
901  return distMin;
902 }
903 
905 //
906 // Calculate distance (<=actual) to closest surface of shape from inside
907 //
909 {
910  G4double rds,roo,roo1, distR, distZ, distMin=0.;
912 
913 #ifdef G4SPECSDEBUG
914  if( Inside(p) == kOutside )
915  {
916  DumpInfo();
917  std::ostringstream message;
918  G4int oldprc = message.precision(16);
919  message << "Point p is outside !?" << G4endl
920  << "Position:" << G4endl
921  << " p.x() = " << p.x()/mm << " mm" << G4endl
922  << " p.y() = " << p.y()/mm << " mm" << G4endl
923  << " p.z() = " << p.z()/mm << " mm";
924  message.precision(oldprc) ;
925  G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
926  JustWarning, message);
927  }
928 #endif
929 
930  // since we have made the above warning, below we are working assuming p
931  // is inside check how close it is to the circular cone with radius equal
932  // to the smaller of the axes
933  //
934  if( sqr(p.x()/minAxis)+sqr(p.y()/minAxis) < sqr(zheight - p.z()) )
935  {
936  rds = std::sqrt(sqr(p.x()) + sqr(p.y()));
937  roo = minAxis*(zheight-p.z()); // radius of cone at z= p.z()
938  roo1 = minAxis*(zheight-zTopCut); // radius of cone at z=+zTopCut
939 
940  distZ=zTopCut - std::fabs(p.z()) ;
941  distR=(roo-rds)/(std::sqrt(1+sqr(minAxis)));
942 
943  if(rds>roo1)
944  {
945  distMin=(zTopCut-p.z())*(roo-rds)/(roo-roo1);
946  distMin=std::min(distMin,distR);
947  }
948  distMin=std::min(distR,distZ);
949  }
950 
951  return distMin;
952 }
953 
955 //
956 // GetEntityType
957 //
959 {
960  return G4String("G4EllipticalCone");
961 }
962 
964 //
965 // Make a clone of the object
966 //
968 {
969  return new G4EllipticalCone(*this);
970 }
971 
973 //
974 // Stream object contents to an output stream
975 //
976 std::ostream& G4EllipticalCone::StreamInfo( std::ostream& os ) const
977 {
978  G4int oldprc = os.precision(16);
979  os << "-----------------------------------------------------------\n"
980  << " *** Dump for solid - " << GetName() << " ***\n"
981  << " ===================================================\n"
982  << " Solid type: G4EllipticalCone\n"
983  << " Parameters: \n"
984 
985  << " semi-axis x: " << xSemiAxis/mm << " mm \n"
986  << " semi-axis y: " << ySemiAxis/mm << " mm \n"
987  << " height z: " << zheight/mm << " mm \n"
988  << " half length in z: " << zTopCut/mm << " mm \n"
989  << "-----------------------------------------------------------\n";
990  os.precision(oldprc);
991 
992  return os;
993 }
994 
996 //
997 // GetPointOnSurface
998 //
999 // returns quasi-uniformly distributed point on surface of elliptical cone
1000 //
1002 {
1003 
1004  G4double phi, sinphi, cosphi, aOne, aTwo, aThree,
1005  chose, zRand, rRand1, rRand2;
1006 
1007  G4double rOne = std::sqrt(sqr(xSemiAxis)
1008  + sqr(ySemiAxis))*(zheight - zTopCut);
1009  G4double rTwo = std::sqrt(sqr(xSemiAxis)
1010  + sqr(ySemiAxis))*(zheight + zTopCut);
1011 
1012  G4int it1=0, it2=0;
1013 
1014  aOne = pi*(rOne + rTwo)*std::sqrt(sqr(rOne - rTwo)+sqr(2.*zTopCut));
1016  aThree = pi*xSemiAxis*ySemiAxis*sqr(zheight-zTopCut);
1017 
1018  phi = RandFlat::shoot(0.,twopi);
1019  cosphi = std::cos(phi);
1020  sinphi = std::sin(phi);
1021 
1022  if(zTopCut >= zheight) aThree = 0.;
1023 
1024  chose = RandFlat::shoot(0.,aOne+aTwo+aThree);
1025  if((chose>=0.) && (chose<aOne))
1026  {
1027  zRand = RandFlat::shoot(-zTopCut,zTopCut);
1028  return G4ThreeVector(xSemiAxis*(zheight-zRand)*cosphi,
1029  ySemiAxis*(zheight-zRand)*sinphi,zRand);
1030  }
1031  else if((chose>=aOne) && (chose<aOne+aTwo))
1032  {
1033  do // Loop checking, 13.08.2015, G.Cosmo
1034  {
1035  rRand1 = RandFlat::shoot(0.,1.) ;
1036  rRand2 = RandFlat::shoot(0.,1.) ;
1037  } while (( rRand2 >= rRand1 ) && (++it1 < 1000)) ;
1038 
1039  // rRand2 = RandFlat::shoot(0.,std::sqrt(1.-sqr(rRand1)));
1040  return G4ThreeVector(rRand1*xSemiAxis*(zheight+zTopCut)*cosphi,
1041  rRand1*ySemiAxis*(zheight+zTopCut)*sinphi, -zTopCut);
1042 
1043  }
1044  // else
1045  //
1046 
1047  do // Loop checking, 13.08.2015, G.Cosmo
1048  {
1049  rRand1 = RandFlat::shoot(0.,1.) ;
1050  rRand2 = RandFlat::shoot(0.,1.) ;
1051  } while (( rRand2 >= rRand1 ) && (++it2 < 1000));
1052 
1053  return G4ThreeVector(rRand1*xSemiAxis*(zheight-zTopCut)*cosphi,
1054  rRand1*ySemiAxis*(zheight-zTopCut)*sinphi, zTopCut);
1055 }
1056 
1057 //
1058 // Methods for visualisation
1059 //
1060 
1062 {
1063  scene.AddSolid(*this);
1064 }
1065 
1067 {
1068  // Define the sides of the box into which the solid instance would fit.
1069  //
1070  G4double maxDim;
1071  maxDim = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
1072  maxDim = maxDim > zTopCut ? maxDim : zTopCut;
1073 
1074  return G4VisExtent (-maxDim, maxDim,
1075  -maxDim, maxDim,
1076  -maxDim, maxDim);
1077 }
1078 
1080 {
1082 }
1083 
1085 {
1086  if ( (!fpPolyhedron)
1089  fpPolyhedron->GetNumberOfRotationSteps()) )
1090  {
1091  G4AutoLock l(&polyhedronMutex);
1092  delete fpPolyhedron;
1094  fRebuildPolyhedron = false;
1095  l.unlock();
1096  }
1097  return fpPolyhedron;
1098 }
G4String GetName() const
G4VSolid * Clone() const
ThreeVector shoot(const G4int Ap, const G4int Af)
static c2_factory< G4double > c2
G4EllipticalCone & operator=(const G4EllipticalCone &rhs)
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4bool GetExtent(G4double &min, G4double &max) const
EInside Inside(const G4ThreeVector &p) const
void SetNormal(const G4ThreeVector &newNormal)
virtual G4bool PartialClip(const G4VoxelLimits &voxelLimit, const EAxis IgnoreMe)
virtual void AddVertexInOrder(const G4ThreeVector vertex)
G4double a
Definition: TRTMaterials.hh:39
virtual void ClearAllVertices()
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
void DumpInfo() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
std::ostream & StreamInfo(std::ostream &os) const
G4ThreeVector GetPointOnSurface() const
double A(double temperature)
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)
static const double twopi
Definition: G4SIunits.hh:75
void AddSurface(const G4ClippablePolygon &surface)
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4int G4Mutex
Definition: G4Threading.hh:173
void SetSemiAxis(G4double x, G4double y, G4double z)
static const double pi
Definition: G4SIunits.hh:74
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
const G4double x[NPOINTSGL]
void SetZCut(G4double newzTopCut)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4Polyhedron * GetPolyhedron() const
G4GeometryType GetEntityType() const
G4VisExtent GetExtent() const
virtual ~G4EllipticalCone()
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
static const double mm
Definition: G4SIunits.hh:114
static G4GeometryTolerance * GetInstance()
void ApplyPointTransform(G4ThreeVector &vec) const
const G4int kMaxMeshSections
Definition: meshdefs.hh:46