Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4EllipticalCone.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/specific/src/G4EllipticalCone.cc (Version 11.3.0) and /geometry/solids/specific/src/G4EllipticalCone.cc (Version 11.0)


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