Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4Ellipsoid.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/G4Ellipsoid.cc (Version 11.3.0) and /geometry/solids/specific/src/G4Ellipsoid.cc (Version 10.2.p3)


  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 // $Id: G4Ellipsoid.cc 102297 2017-01-20 13:33:54Z gcosmo $
                                                   >>  27 //
 26 // class G4Ellipsoid                               28 // class G4Ellipsoid
 27 //                                                 29 //
 28 // Implementation of G4Ellipsoid class         <<  30 // Implementation for G4Ellipsoid class
                                                   >>  31 //
                                                   >>  32 // History:
                                                   >>  33 //
                                                   >>  34 // 10.11.99 G.Horton-Smith  -- first writing, based on G4Sphere class
                                                   >>  35 // 25.02.05 G.Guerrieri -- Modified for future Geant4 release
 29 //                                                 36 //
 30 // 10.11.99 G.Horton-Smith: first writing, bas << 
 31 // 25.02.05 G.Guerrieri: Revised               << 
 32 // 15.12.19 E.Tcherniaev: Complete revision    << 
 33 // -------------------------------------------     37 // --------------------------------------------------------------------
 34                                                    38 
 35 #include "G4Ellipsoid.hh"                      << 
 36                                                << 
 37 #if !(defined(G4GEOM_USE_UELLIPSOID) && define << 
 38                                                << 
 39 #include "globals.hh"                              39 #include "globals.hh"
 40                                                    40 
                                                   >>  41 #include "G4Ellipsoid.hh"
                                                   >>  42 
 41 #include "G4VoxelLimits.hh"                        43 #include "G4VoxelLimits.hh"
 42 #include "G4AffineTransform.hh"                    44 #include "G4AffineTransform.hh"
 43 #include "G4GeometryTolerance.hh"                  45 #include "G4GeometryTolerance.hh"
 44 #include "G4BoundingEnvelope.hh"               <<  46 
 45 #include "G4RandomTools.hh"                    <<  47 #include "meshdefs.hh"
 46 #include "G4QuickRand.hh"                      <<  48 #include "Randomize.hh"
 47                                                    49 
 48 #include "G4VPVParameterisation.hh"                50 #include "G4VPVParameterisation.hh"
 49                                                    51 
 50 #include "G4VGraphicsScene.hh"                     52 #include "G4VGraphicsScene.hh"
 51 #include "G4VisExtent.hh"                          53 #include "G4VisExtent.hh"
 52                                                    54 
 53 #include "G4AutoLock.hh"                           55 #include "G4AutoLock.hh"
 54                                                    56 
 55 namespace                                          57 namespace
 56 {                                                  58 {
 57   G4Mutex polyhedronMutex  = G4MUTEX_INITIALIZ <<  59   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
 58   G4Mutex lateralareaMutex = G4MUTEX_INITIALIZ << 
 59 }                                                  60 }
 60                                                    61 
 61 using namespace CLHEP;                             62 using namespace CLHEP;
 62                                                    63 
 63 ////////////////////////////////////////////// <<  64 ///////////////////////////////////////////////////////////////////////////////
 64 //                                                 65 //
 65 // Constructor                                 <<  66 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
                                                   >>  67 //             - note if pDPhi>2PI then reset to 2PI
 66                                                    68 
 67 G4Ellipsoid::G4Ellipsoid(const G4String& name, <<  69 G4Ellipsoid::G4Ellipsoid(const G4String& pName,
 68                                G4double xSemiA <<  70                                G4double pxSemiAxis,
 69                                G4double ySemiA <<  71                                G4double pySemiAxis,
 70                                G4double zSemiA <<  72                                G4double pzSemiAxis,
 71                                G4double zBotto <<  73                                G4double pzBottomCut,
 72                                G4double zTopCu <<  74                                G4double pzTopCut)
 73   : G4VSolid(name), fDx(xSemiAxis), fDy(ySemiA <<  75   : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
 74     fZBottomCut(zBottomCut), fZTopCut(zTopCut) <<  76     fCubicVolume(0.), fSurfaceArea(0.), zBottomCut(0.), zTopCut(0.)
 75 {                                              <<  77 {
 76   CheckParameters();                           <<  78   // note: for users that want to use the full ellipsoid it is useful
                                                   >>  79   // to include a default for the cuts 
                                                   >>  80 
                                                   >>  81   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
                                                   >>  82 
                                                   >>  83   halfCarTolerance = kCarTolerance*0.5;
                                                   >>  84   halfRadTolerance = kRadTolerance*0.5;
                                                   >>  85 
                                                   >>  86   // Check Semi-Axis
                                                   >>  87   if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
                                                   >>  88   {
                                                   >>  89      std::ostringstream message;
                                                   >>  90      message << "Invalid semi-axis - " << GetName();
                                                   >>  91      G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
                                                   >>  92                  FatalErrorInArgument, message);
                                                   >>  93   }
                                                   >>  94   SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis);
                                                   >>  95 
                                                   >>  96   if ( pzBottomCut == 0 && pzTopCut == 0 )
                                                   >>  97   {
                                                   >>  98      SetZCuts(-pzSemiAxis, pzSemiAxis);
                                                   >>  99   }
                                                   >> 100   else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
                                                   >> 101          && (pzBottomCut < pzTopCut) )
                                                   >> 102   {
                                                   >> 103      SetZCuts(pzBottomCut, pzTopCut);
                                                   >> 104   }
                                                   >> 105   else
                                                   >> 106   {
                                                   >> 107      std::ostringstream message;
                                                   >> 108      message << "Invalid z-coordinate for cutting plane - " << GetName();
                                                   >> 109      G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
                                                   >> 110                  FatalErrorInArgument, message);
                                                   >> 111   }
 77 }                                                 112 }
 78                                                   113 
 79 ////////////////////////////////////////////// << 114 ///////////////////////////////////////////////////////////////////////////////
 80 //                                                115 //
 81 // Fake default constructor - sets only member    116 // Fake default constructor - sets only member data and allocates memory
 82 //                            for usage restri << 117 //                            for usage restricted to object persistency.
 83                                                << 118 //
 84 G4Ellipsoid::G4Ellipsoid( __void__& a )           119 G4Ellipsoid::G4Ellipsoid( __void__& a )
 85   : G4VSolid(a), fDx(0.), fDy(0.), fDz(0.), fZ << 120   : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0), kRadTolerance(0.),
                                                   >> 121     halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
                                                   >> 122     fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
                                                   >> 123     semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
 86 {                                                 124 {
 87 }                                                 125 }
 88                                                   126 
 89 ////////////////////////////////////////////// << 127 ///////////////////////////////////////////////////////////////////////////////
 90 //                                                128 //
 91 // Destructor                                     129 // Destructor
 92                                                   130 
 93 G4Ellipsoid::~G4Ellipsoid()                       131 G4Ellipsoid::~G4Ellipsoid()
 94 {                                                 132 {
 95   delete fpPolyhedron; fpPolyhedron = nullptr; << 133   delete fpPolyhedron; fpPolyhedron = 0;
 96 }                                                 134 }
 97                                                   135 
 98 ////////////////////////////////////////////// << 136 ///////////////////////////////////////////////////////////////////////////////
 99 //                                                137 //
100 // Copy constructor                               138 // Copy constructor
101                                                   139 
102 G4Ellipsoid::G4Ellipsoid(const G4Ellipsoid& rh    140 G4Ellipsoid::G4Ellipsoid(const G4Ellipsoid& rhs)
103   : G4VSolid(rhs),                                141   : G4VSolid(rhs),
104    fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),   << 142     fRebuildPolyhedron(false), fpPolyhedron(0),
105    fZBottomCut(rhs.fZBottomCut), fZTopCut(rhs. << 143     kRadTolerance(rhs.kRadTolerance),
106    halfTolerance(rhs.halfTolerance),           << 144     halfCarTolerance(rhs.halfCarTolerance),
107    fXmax(rhs.fXmax), fYmax(rhs.fYmax),         << 145     halfRadTolerance(rhs.halfRadTolerance),
108    fRsph(rhs.fRsph), fR(rhs.fR),               << 146     fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
109    fSx(rhs.fSx), fSy(rhs.fSy), fSz(rhs.fSz),   << 147     xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
110    fZMidCut(rhs.fZMidCut), fZDimCut(rhs.fZDimC << 148     zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
111    fQ1(rhs.fQ1), fQ2(rhs.fQ2),                 << 149     zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
112    fCubicVolume(rhs.fCubicVolume),             << 
113    fSurfaceArea(rhs.fSurfaceArea),             << 
114    fLateralArea(rhs.fLateralArea)              << 
115 {                                                 150 {
116 }                                                 151 }
117                                                   152 
118 ////////////////////////////////////////////// << 153 ///////////////////////////////////////////////////////////////////////////////
119 //                                                154 //
120 // Assignment operator                            155 // Assignment operator
121                                                   156 
122 G4Ellipsoid& G4Ellipsoid::operator = (const G4 << 157 G4Ellipsoid& G4Ellipsoid::operator = (const G4Ellipsoid& rhs) 
123 {                                                 158 {
124    // Check assignment to self                    159    // Check assignment to self
125    //                                             160    //
126    if (this == &rhs)  { return *this; }           161    if (this == &rhs)  { return *this; }
127                                                   162 
128    // Copy base class data                        163    // Copy base class data
129    //                                             164    //
130    G4VSolid::operator=(rhs);                      165    G4VSolid::operator=(rhs);
131                                                   166 
132    // Copy data                                   167    // Copy data
133    //                                             168    //
134    fDx = rhs.fDx;                              << 169    kRadTolerance = rhs.kRadTolerance;
135    fDy = rhs.fDy;                              << 170    halfCarTolerance = rhs.halfCarTolerance;
136    fDz = rhs.fDz;                              << 171    halfRadTolerance = rhs.halfRadTolerance;
137    fZBottomCut = rhs.fZBottomCut;              << 172    fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
138    fZTopCut = rhs.fZTopCut;                    << 173    xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
139                                                << 174    zSemiAxis = rhs.zSemiAxis; semiAxisMax = rhs.semiAxisMax;
140    halfTolerance = rhs.halfTolerance;          << 175    zBottomCut = rhs.zBottomCut; zTopCut = rhs.zTopCut;
141    fXmax = rhs.fXmax;                          << 
142    fYmax = rhs.fYmax;                          << 
143    fRsph = rhs.fRsph;                          << 
144    fR = rhs.fR;                                << 
145    fSx = rhs.fSx;                              << 
146    fSy = rhs.fSy;                              << 
147    fSz = rhs.fSz;                              << 
148    fZMidCut = rhs.fZMidCut;                    << 
149    fZDimCut = rhs.fZDimCut;                    << 
150    fQ1 = rhs.fQ1;                              << 
151    fQ2 = rhs.fQ2;                              << 
152                                                << 
153    fCubicVolume = rhs.fCubicVolume;            << 
154    fSurfaceArea = rhs.fSurfaceArea;            << 
155    fLateralArea = rhs.fLateralArea;            << 
156                                                << 
157    fRebuildPolyhedron = false;                    176    fRebuildPolyhedron = false;
158    delete fpPolyhedron; fpPolyhedron = nullptr << 177    delete fpPolyhedron; fpPolyhedron = 0;
159                                                   178 
160    return *this;                                  179    return *this;
161 }                                                 180 }
162                                                   181 
163 ////////////////////////////////////////////// << 182 ////////////////////////////////////////////////////////////////////////
164 //                                             << 
165 // Check parameters and make precalculation    << 
166                                                << 
167 void G4Ellipsoid::CheckParameters()            << 
168 {                                              << 
169   halfTolerance = 0.5 * kCarTolerance; // half << 
170   G4double dmin = 2. * kCarTolerance;          << 
171                                                << 
172   // Check dimensions                          << 
173   //                                           << 
174   if (fDx < dmin || fDy < dmin || fDz < dmin)  << 
175   {                                            << 
176     std::ostringstream message;                << 
177     message << "Invalid (too small or negative << 
178             << GetName()  << "\n"              << 
179             << "  semi-axis x: " << fDx << "\n << 
180             << "  semi-axis y: " << fDy << "\n << 
181             << "  semi-axis z: " << fDz;       << 
182     G4Exception("G4Ellipsoid::CheckParameters( << 
183                 FatalException, message);      << 
184   }                                            << 
185   G4double A = fDx;                            << 
186   G4double B = fDy;                            << 
187   G4double C = fDz;                            << 
188                                                << 
189   // Check cuts                                << 
190   //                                           << 
191   if (fZBottomCut == 0. && fZTopCut == 0.)     << 
192   {                                            << 
193     fZBottomCut = -C;                          << 
194     fZTopCut = C;                              << 
195   }                                            << 
196   if (fZBottomCut >= C || fZTopCut <= -C || fZ << 
197   {                                            << 
198     std::ostringstream message;                << 
199     message << "Invalid Z cuts for Solid: "    << 
200             << GetName() << "\n"               << 
201             << "  bottom cut: " << fZBottomCut << 
202             << "  top cut: " << fZTopCut;      << 
203     G4Exception("G4Ellipsoid::CheckParameters( << 
204                 FatalException, message);      << 
205                                                << 
206   }                                            << 
207   fZBottomCut = std::max(fZBottomCut, -C);     << 
208   fZTopCut = std::min(fZTopCut, C);            << 
209                                                << 
210   // Set extent in x and y                     << 
211   fXmax = A;                                   << 
212   fYmax = B;                                   << 
213   if (fZBottomCut > 0.)                        << 
214   {                                            << 
215     G4double ratio = fZBottomCut / C;          << 
216     G4double scale = std::sqrt((1. - ratio) *  << 
217     fXmax *= scale;                            << 
218     fYmax *= scale;                            << 
219   }                                            << 
220   if (fZTopCut < 0.)                           << 
221   {                                            << 
222     G4double ratio  = fZTopCut / C;            << 
223     G4double scale  = std::sqrt((1. - ratio) * << 
224     fXmax *= scale;                            << 
225     fYmax *= scale;                            << 
226   }                                            << 
227                                                << 
228   // Set scale factors                         << 
229   fRsph = std::max(std::max(A, B), C); // boun << 
230   fR    = std::min(std::min(A, B), C); // radi << 
231   fSx   = fR / A; // X scale factor            << 
232   fSy   = fR / B; // Y scale factor            << 
233   fSz   = fR / C; // Z scale factor            << 
234                                                << 
235   // Scaled cuts                               << 
236   fZMidCut = 0.5 * (fZTopCut + fZBottomCut) *  << 
237   fZDimCut = 0.5 * (fZTopCut - fZBottomCut) *  << 
238                                                << 
239   // Coefficients for approximation of distanc << 
240   fQ1 = 0.5 / fR;                              << 
241   fQ2 = 0.5 * fR + halfTolerance * halfToleran << 
242                                                << 
243   fCubicVolume = 0.; // volume                 << 
244   fSurfaceArea = 0.; // surface area           << 
245   fLateralArea = 0.; // lateral surface area   << 
246 }                                              << 
247                                                << 
248 ////////////////////////////////////////////// << 
249 //                                                183 //
250 // Dispatch to parameterisation for replicatio    184 // Dispatch to parameterisation for replication mechanism dimension
251 // computation & modification                  << 185 // computation & modification.
252                                                   186 
253 void G4Ellipsoid::ComputeDimensions(G4VPVParam    187 void G4Ellipsoid::ComputeDimensions(G4VPVParameterisation* p,
254                                     const G4in    188                                     const G4int n,
255                                     const G4VP    189                                     const G4VPhysicalVolume* pRep)
256 {                                                 190 {
257   p->ComputeDimensions(*this,n,pRep);             191   p->ComputeDimensions(*this,n,pRep);
258 }                                                 192 }
259                                                   193 
260 ////////////////////////////////////////////// << 194 ///////////////////////////////////////////////////////////////////////////////
261 //                                             << 
262 // Get bounding box                            << 
263                                                << 
264 void G4Ellipsoid::BoundingLimits(G4ThreeVector << 
265                                  G4ThreeVector << 
266 {                                              << 
267   pMin.set(-fXmax,-fYmax, fZBottomCut);        << 
268   pMax.set( fXmax, fYmax, fZTopCut);           << 
269 }                                              << 
270                                                << 
271 ////////////////////////////////////////////// << 
272 //                                                195 //
273 // Calculate extent under transform and specif << 196 // Calculate extent under transform and specified limit
274                                                   197 
275 G4bool                                            198 G4bool
276 G4Ellipsoid::CalculateExtent(const EAxis pAxis    199 G4Ellipsoid::CalculateExtent(const EAxis pAxis,
277                              const G4VoxelLimi    200                              const G4VoxelLimits& pVoxelLimit,
278                              const G4AffineTra    201                              const G4AffineTransform& pTransform,
279                                    G4double& p    202                                    G4double& pMin, G4double& pMax) const
280 {                                                 203 {
281   G4ThreeVector bmin, bmax;                    << 204   if (!pTransform.IsRotated())
                                                   >> 205   {
                                                   >> 206     // Special case handling for unrotated solid ellipsoid
                                                   >> 207     // Compute x/y/z mins and maxs for bounding box respecting limits,
                                                   >> 208     // with early returns if outside limits. Then switch() on pAxis,
                                                   >> 209     // and compute exact x and y limit for x/y case
                                                   >> 210 
                                                   >> 211     G4double xoffset,xMin,xMax;
                                                   >> 212     G4double yoffset,yMin,yMax;
                                                   >> 213     G4double zoffset,zMin,zMax;
                                                   >> 214 
                                                   >> 215     G4double maxDiff,newMin,newMax;
                                                   >> 216     G4double xoff,yoff;
                                                   >> 217 
                                                   >> 218     xoffset=pTransform.NetTranslation().x();
                                                   >> 219     xMin=xoffset - xSemiAxis;
                                                   >> 220     xMax=xoffset + xSemiAxis;
                                                   >> 221     if (pVoxelLimit.IsXLimited())
                                                   >> 222     {
                                                   >> 223       if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
                                                   >> 224         || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
                                                   >> 225       {
                                                   >> 226         return false;
                                                   >> 227       }
                                                   >> 228       else
                                                   >> 229       {
                                                   >> 230         if (xMin<pVoxelLimit.GetMinXExtent())
                                                   >> 231         {
                                                   >> 232           xMin=pVoxelLimit.GetMinXExtent();
                                                   >> 233         }
                                                   >> 234         if (xMax>pVoxelLimit.GetMaxXExtent())
                                                   >> 235         {
                                                   >> 236           xMax=pVoxelLimit.GetMaxXExtent();
                                                   >> 237         }
                                                   >> 238       }
                                                   >> 239     }
                                                   >> 240 
                                                   >> 241     yoffset=pTransform.NetTranslation().y();
                                                   >> 242     yMin=yoffset - ySemiAxis;
                                                   >> 243     yMax=yoffset + ySemiAxis;
                                                   >> 244     if (pVoxelLimit.IsYLimited())
                                                   >> 245     {
                                                   >> 246       if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
                                                   >> 247         || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
                                                   >> 248       {
                                                   >> 249         return false;
                                                   >> 250       }
                                                   >> 251       else
                                                   >> 252       {
                                                   >> 253         if (yMin<pVoxelLimit.GetMinYExtent())
                                                   >> 254         {
                                                   >> 255           yMin=pVoxelLimit.GetMinYExtent();
                                                   >> 256         }
                                                   >> 257         if (yMax>pVoxelLimit.GetMaxYExtent())
                                                   >> 258         {
                                                   >> 259           yMax=pVoxelLimit.GetMaxYExtent();
                                                   >> 260         }
                                                   >> 261       }
                                                   >> 262     }
282                                                   263 
283   // Get bounding box                          << 264     zoffset=pTransform.NetTranslation().z();
284   BoundingLimits(bmin,bmax);                   << 265     zMin=zoffset + (-zSemiAxis > zBottomCut ? -zSemiAxis : zBottomCut);
                                                   >> 266     zMax=zoffset + ( zSemiAxis < zTopCut ? zSemiAxis : zTopCut);
                                                   >> 267     if (pVoxelLimit.IsZLimited())
                                                   >> 268     {
                                                   >> 269       if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
                                                   >> 270         || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
                                                   >> 271       {
                                                   >> 272         return false;
                                                   >> 273       }
                                                   >> 274       else
                                                   >> 275       {
                                                   >> 276         if (zMin<pVoxelLimit.GetMinZExtent())
                                                   >> 277         {
                                                   >> 278           zMin=pVoxelLimit.GetMinZExtent();
                                                   >> 279         }
                                                   >> 280         if (zMax>pVoxelLimit.GetMaxZExtent())
                                                   >> 281         {
                                                   >> 282           zMax=pVoxelLimit.GetMaxZExtent();
                                                   >> 283         }
                                                   >> 284       }
                                                   >> 285     }
285                                                   286 
286   // Find extent                               << 287     // if here, then known to cut bounding box around ellipsoid
287   G4BoundingEnvelope bbox(bmin,bmax);          << 288     //
288   return bbox.CalculateExtent(pAxis,pVoxelLimi << 289     xoff = (xoffset < xMin) ? (xMin-xoffset)
289 }                                              << 290          : (xoffset > xMax) ? (xoffset-xMax) : 0.0;
                                                   >> 291     yoff = (yoffset < yMin) ? (yMin-yoffset)
                                                   >> 292          : (yoffset > yMax) ? (yoffset-yMax) : 0.0;
                                                   >> 293 
                                                   >> 294     // detailed calculations
                                                   >> 295     // NOTE: does not use X or Y offsets to adjust Z range,
                                                   >> 296     // and does not use Z offset to adjust X or Y range,
                                                   >> 297     // which is consistent with G4Sphere::CalculateExtent behavior
                                                   >> 298     //
                                                   >> 299     switch (pAxis)
                                                   >> 300     {
                                                   >> 301       case kXAxis:
                                                   >> 302         if (yoff==0.)
                                                   >> 303         {
                                                   >> 304           // YZ limits cross max/min x => no change
                                                   >> 305           //
                                                   >> 306           pMin=xMin;
                                                   >> 307           pMax=xMax;
                                                   >> 308         }
                                                   >> 309         else
                                                   >> 310         {
                                                   >> 311           // YZ limits don't cross max/min x => compute max delta x,
                                                   >> 312           // hence new mins/maxs
                                                   >> 313           //
                                                   >> 314           maxDiff= 1.0-sqr(yoff/ySemiAxis);
                                                   >> 315           if (maxDiff < 0.0) { return false; }
                                                   >> 316           maxDiff= xSemiAxis * std::sqrt(maxDiff);
                                                   >> 317           newMin=xoffset-maxDiff;
                                                   >> 318           newMax=xoffset+maxDiff;
                                                   >> 319           pMin=(newMin<xMin) ? xMin : newMin;
                                                   >> 320           pMax=(newMax>xMax) ? xMax : newMax;
                                                   >> 321         }
                                                   >> 322         break;
                                                   >> 323       case kYAxis:
                                                   >> 324         if (xoff==0.)
                                                   >> 325         {
                                                   >> 326           // XZ limits cross max/min y => no change
                                                   >> 327           //
                                                   >> 328           pMin=yMin;
                                                   >> 329           pMax=yMax;
                                                   >> 330         }
                                                   >> 331         else
                                                   >> 332         {
                                                   >> 333           // XZ limits don't cross max/min y => compute max delta y,
                                                   >> 334           // hence new mins/maxs
                                                   >> 335           //
                                                   >> 336           maxDiff= 1.0-sqr(xoff/xSemiAxis);
                                                   >> 337           if (maxDiff < 0.0) { return false; }
                                                   >> 338           maxDiff= ySemiAxis * std::sqrt(maxDiff);
                                                   >> 339           newMin=yoffset-maxDiff;
                                                   >> 340           newMax=yoffset+maxDiff;
                                                   >> 341           pMin=(newMin<yMin) ? yMin : newMin;
                                                   >> 342           pMax=(newMax>yMax) ? yMax : newMax;
                                                   >> 343         }
                                                   >> 344         break;
                                                   >> 345       case kZAxis:
                                                   >> 346         pMin=zMin;
                                                   >> 347         pMax=zMax;
                                                   >> 348         break;
                                                   >> 349       default:
                                                   >> 350         break;
                                                   >> 351     }
                                                   >> 352   
                                                   >> 353     pMin-=kCarTolerance;
                                                   >> 354     pMax+=kCarTolerance;
                                                   >> 355     return true;
                                                   >> 356   }
                                                   >> 357   else  // not rotated
                                                   >> 358   {
                                                   >> 359     G4int i,j,noEntries,noBetweenSections;
                                                   >> 360     G4bool existsAfterClip=false;
290                                                   361 
291 ////////////////////////////////////////////// << 362     // Calculate rotated vertex coordinates
292 //                                             << 
293 // Return position of point: inside/outside/on << 
294                                                   363 
295 EInside G4Ellipsoid::Inside(const G4ThreeVecto << 364     G4int noPolygonVertices=0;
296 {                                              << 365     G4ThreeVectorList* vertices =
297   G4double x     = p.x() * fSx;                << 366       CreateRotatedVertices(pTransform,noPolygonVertices);
298   G4double y     = p.y() * fSy;                << 367 
299   G4double z     = p.z() * fSz;                << 368     pMin=+kInfinity;
300   G4double rr    = x * x + y * y + z * z;      << 369     pMax=-kInfinity;
301   G4double distZ = std::abs(z - fZMidCut) - fZ << 370 
302   G4double distR = fQ1 * rr - fQ2;             << 371     noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections
303   G4double dist  = std::max(distZ, distR);     << 372     noBetweenSections=noEntries-noPolygonVertices;
                                                   >> 373     
                                                   >> 374     G4ThreeVectorList ThetaPolygon;
                                                   >> 375     for (i=0;i<noEntries;i+=noPolygonVertices)
                                                   >> 376     {
                                                   >> 377       for(j=0;j<(noPolygonVertices/2)-1;j++)
                                                   >> 378       {
                                                   >> 379         ThetaPolygon.push_back((*vertices)[i+j]);  
                                                   >> 380         ThetaPolygon.push_back((*vertices)[i+j+1]);  
                                                   >> 381         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]);
                                                   >> 382         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]);
                                                   >> 383         CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 384         ThetaPolygon.clear();
                                                   >> 385       }
                                                   >> 386     }
                                                   >> 387     for (i=0;i<noBetweenSections;i+=noPolygonVertices)
                                                   >> 388     {
                                                   >> 389       for(j=0;j<noPolygonVertices-1;j++)
                                                   >> 390       {
                                                   >> 391         ThetaPolygon.push_back((*vertices)[i+j]);  
                                                   >> 392         ThetaPolygon.push_back((*vertices)[i+j+1]);  
                                                   >> 393         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]);
                                                   >> 394         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]);
                                                   >> 395         CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 396         ThetaPolygon.clear();
                                                   >> 397       }
                                                   >> 398       ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]);
                                                   >> 399       ThetaPolygon.push_back((*vertices)[i]);
                                                   >> 400       ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]);
                                                   >> 401       ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]);
                                                   >> 402       CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 403       ThetaPolygon.clear();
                                                   >> 404     }
                                                   >> 405     if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
                                                   >> 406     {
                                                   >> 407       existsAfterClip=true;
                                                   >> 408     
                                                   >> 409       // Add 2*tolerance to avoid precision troubles
                                                   >> 410       //
                                                   >> 411       pMin-=kCarTolerance;
                                                   >> 412       pMax+=kCarTolerance;
304                                                   413 
305   if (dist > halfTolerance) return kOutside;   << 414     }
306   return (dist > -halfTolerance) ? kSurface :  << 415     else
                                                   >> 416     {
                                                   >> 417       // Check for case where completely enveloping clipping volume
                                                   >> 418       // If point inside then we are confident that the solid completely
                                                   >> 419       // envelopes the clipping volume. Hence set min/max extents according
                                                   >> 420       // to clipping volume extents along the specified axis.
                                                   >> 421       //
                                                   >> 422       G4ThreeVector
                                                   >> 423       clipCentre((pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
                                                   >> 424                  (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
                                                   >> 425                  (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
                                                   >> 426 
                                                   >> 427       if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
                                                   >> 428       {
                                                   >> 429         existsAfterClip=true;
                                                   >> 430         pMin=pVoxelLimit.GetMinExtent(pAxis);
                                                   >> 431         pMax=pVoxelLimit.GetMaxExtent(pAxis);
                                                   >> 432       }
                                                   >> 433     }
                                                   >> 434     delete vertices;
                                                   >> 435     return existsAfterClip;
                                                   >> 436   }
307 }                                                 437 }
308                                                   438 
309 ////////////////////////////////////////////// << 439 ///////////////////////////////////////////////////////////////////////////////
310 //                                                440 //
311 // Return unit normal to surface at p          << 441 // Return whether point inside/outside/on surface
                                                   >> 442 // Split into radius, phi, theta checks
                                                   >> 443 // Each check modifies `in', or returns as approprate
312                                                   444 
313 G4ThreeVector G4Ellipsoid::SurfaceNormal( cons << 445 EInside G4Ellipsoid::Inside(const G4ThreeVector& p) const
314 {                                                 446 {
315   G4ThreeVector norm(0., 0., 0.);              << 447   G4double rad2oo,  // outside surface outer tolerance
316   G4int nsurf = 0;                             << 448            rad2oi;  // outside surface inner tolerance
                                                   >> 449   EInside in;
317                                                   450 
318   // Check cuts                                << 451   // check this side of z cut first, because that's fast
319   G4double x = p.x() * fSx;                    << 452   //
320   G4double y = p.y() * fSy;                    << 453   if (p.z() < zBottomCut-halfRadTolerance) { return in=kOutside; }
321   G4double z = p.z() * fSz;                    << 454   if (p.z() > zTopCut+halfRadTolerance)    { return in=kOutside; }
322   G4double distZ = std::abs(z - fZMidCut) - fZ << 
323   if (std::abs(distZ) <= halfTolerance)        << 
324   {                                            << 
325     norm.setZ(std::copysign(1., z - fZMidCut)) << 
326     ++nsurf;                                   << 
327   }                                            << 
328                                                   455 
329   // Check lateral surface                     << 456   rad2oo= sqr(p.x()/(xSemiAxis+halfRadTolerance))
330   G4double distR = fQ1*(x*x + y*y + z*z) - fQ2 << 457         + sqr(p.y()/(ySemiAxis+halfRadTolerance))
331   if (std::abs(distR) <= halfTolerance)        << 458         + sqr(p.z()/(zSemiAxis+halfRadTolerance));
                                                   >> 459 
                                                   >> 460   if (rad2oo > 1.0)  { return in=kOutside; }
                                                   >> 461     
                                                   >> 462   rad2oi= sqr(p.x()*(1.0+halfRadTolerance/xSemiAxis)/xSemiAxis)
                                                   >> 463       + sqr(p.y()*(1.0+halfRadTolerance/ySemiAxis)/ySemiAxis)
                                                   >> 464       + sqr(p.z()*(1.0+halfRadTolerance/zSemiAxis)/zSemiAxis);
                                                   >> 465 
                                                   >> 466   // Check radial surfaces
                                                   >> 467   //  sets `in' (already checked for rad2oo > 1.0)
                                                   >> 468   //
                                                   >> 469   if (rad2oi < 1.0)
332   {                                               470   {
333     // normal = (p.x/a^2, p.y/b^2, p.z/c^2)    << 471     in = ( (p.z() < zBottomCut+halfRadTolerance)
334     norm += G4ThreeVector(x*fSx, y*fSy, z*fSz) << 472         || (p.z() > zTopCut-halfRadTolerance) ) ? kSurface : kInside;
335     ++nsurf;                                   << 473     if ( rad2oi > 1.0-halfRadTolerance )  { in=kSurface; }
336   }                                               474   }
337                                                << 475   else 
338   // Return normal                             << 
339   if (nsurf == 1) return norm;                 << 
340   else if (nsurf > 1) return norm.unit(); // e << 
341   else                                         << 
342   {                                               476   {
343 #ifdef G4SPECSDEBUG                            << 477     in = kSurface;
344     std::ostringstream message;                << 
345     G4long oldprc = message.precision(16);     << 
346     message << "Point p is not on surface (!?) << 
347             << GetName() << "\n";              << 
348     message << "Position:\n";                  << 
349     message << "   p.x() = " << p.x()/mm << "  << 
350     message << "   p.y() = " << p.y()/mm << "  << 
351     message << "   p.z() = " << p.z()/mm << "  << 
352     G4cout.precision(oldprc);                  << 
353     G4Exception("G4Ellipsoid::SurfaceNormal(p) << 
354                 JustWarning, message );        << 
355     DumpInfo();                                << 
356 #endif                                         << 
357     return ApproxSurfaceNormal(p);             << 
358   }                                               478   }
                                                   >> 479   return in;
                                                   >> 480 
359 }                                                 481 }
360                                                   482 
361 ////////////////////////////////////////////// << 483 ///////////////////////////////////////////////////////////////////////////////
362 //                                                484 //
363 // Find surface nearest to point and return co << 485 // Return unit normal of surface closest to p not protected against p=0
364 // This method normally should not be called.  << 
365                                                   486 
366 G4ThreeVector G4Ellipsoid::ApproxSurfaceNormal << 487 G4ThreeVector G4Ellipsoid::SurfaceNormal( const G4ThreeVector& p) const
367 {                                                 488 {
368   G4double x  = p.x() * fSx;                   << 489   G4double distR, distZBottom, distZTop;
369   G4double y  = p.y() * fSy;                   << 490 
370   G4double z  = p.z() * fSz;                   << 491   // normal vector with special magnitude:  parallel to normal, units 1/length
371   G4double rr = x * x + y * y + z * z;         << 492   // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside
372   G4double distZ = std::abs(z - fZMidCut) - fZ << 493   //
373   G4double distR = std::sqrt(rr) - fR;         << 494   G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
374   if  (distR > distZ && rr > 0.) // distR > di << 495                      p.y()/(ySemiAxis*ySemiAxis),
375     return G4ThreeVector(x*fSx, y*fSy, z*fSz). << 496                      p.z()/(zSemiAxis*zSemiAxis));
376   else                                         << 497   G4double radius = 1.0/norm.mag();
377     return { 0., 0., std::copysign(1., z - fZM << 498 
                                                   >> 499   // approximate distance to curved surface
                                                   >> 500   //
                                                   >> 501   distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
                                                   >> 502 
                                                   >> 503   // Distance to z-cut plane
                                                   >> 504   //
                                                   >> 505   distZBottom = std::fabs( p.z() - zBottomCut );
                                                   >> 506   distZTop = std::fabs( p.z() - zTopCut );
                                                   >> 507 
                                                   >> 508   if ( (distZBottom < distR) || (distZTop < distR) )
                                                   >> 509   {
                                                   >> 510     return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
                                                   >> 511   }
                                                   >> 512   return ( norm *= radius );
378 }                                                 513 }
379                                                   514 
380 ////////////////////////////////////////////// << 515 ///////////////////////////////////////////////////////////////////////////////
                                                   >> 516 //
                                                   >> 517 // Calculate distance to shape from outside, along normalised vector
                                                   >> 518 // - return kInfinity if no intersection, or intersection distance <= tolerance
381 //                                                519 //
382 // Calculate distance to shape from outside al << 
383                                                   520 
384 G4double G4Ellipsoid::DistanceToIn(const G4Thr << 521 G4double G4Ellipsoid::DistanceToIn( const G4ThreeVector& p,
385                                    const G4Thr << 522                                     const G4ThreeVector& v  ) const
386 {                                                 523 {
387   G4double offset = 0.;                        << 524   G4double distMin = std::min(xSemiAxis,ySemiAxis);
388   G4ThreeVector pcur = p;                      << 525   const G4double dRmax = 100.*std::min(distMin,zSemiAxis);
389                                                << 526   distMin= kInfinity;
390   // Check if point is flying away, relative t << 
391   //                                           << 
392   G4double safex = std::abs(p.x()) - fXmax;    << 
393   G4double safey = std::abs(p.y()) - fYmax;    << 
394   G4double safet = p.z() - fZTopCut;           << 
395   G4double safeb = fZBottomCut - p.z();        << 
396                                                   527 
397   if (safex >= -halfTolerance && p.x() * v.x() << 528   // check to see if Z plane is relevant
398   if (safey >= -halfTolerance && p.y() * v.y() << 529   if (p.z() <= zBottomCut+halfCarTolerance)
399   if (safet >= -halfTolerance && v.z() >= 0.)  << 530   {
400   if (safeb >= -halfTolerance && v.z() <= 0.)  << 531     if (v.z() <= 0.0) { return distMin; }
                                                   >> 532     G4double distZ = (zBottomCut - p.z()) / v.z();
401                                                   533 
402   // Relocate point, if required               << 534     if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
403   //                                           << 535     {
404   G4double safe = std::max(std::max(std::max(s << 536       // early exit since can't intercept curved surface if we reach here
405   if (safe > 32. * fRsph)                      << 537       if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
                                                   >> 538       return distMin= distZ;
                                                   >> 539     }
                                                   >> 540   }
                                                   >> 541   if (p.z() >= zTopCut-halfCarTolerance)
406   {                                               542   {
407     offset = (1. - 1.e-08) * safe - 2. * fRsph << 543     if (v.z() >= 0.0) { return distMin;}
408     pcur += offset * v;                        << 544     G4double distZ = (zTopCut - p.z()) / v.z();
409     G4double dist = DistanceToIn(pcur, v);     << 545     if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
410     return (dist == kInfinity) ? kInfinity : d << 546     {
                                                   >> 547       // early exit since can't intercept curved surface if we reach here
                                                   >> 548       if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
                                                   >> 549       return distMin= distZ;
                                                   >> 550     }
411   }                                               551   }
                                                   >> 552   // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface
412                                                   553 
413   // Scale ellipsoid to sphere                 << 554   // now check curved surface intercept
414   //                                           << 555   G4double A,B,C;
415   G4double px = pcur.x() * fSx;                << 
416   G4double py = pcur.y() * fSy;                << 
417   G4double pz = pcur.z() * fSz;                << 
418   G4double vx = v.x() * fSx;                   << 
419   G4double vy = v.y() * fSy;                   << 
420   G4double vz = v.z() * fSz;                   << 
421                                                   556 
422   // Check if point is leaving the solid       << 557   A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
423   //                                           << 558   C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0;
424   G4double dzcut = fZDimCut;                   << 559   B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis)
425   G4double pzcut = pz - fZMidCut;              << 560            + p.y()*v.y()/(ySemiAxis*ySemiAxis)
426   G4double distZ = std::abs(pzcut) - dzcut;    << 561            + p.z()*v.z()/(zSemiAxis*zSemiAxis) );
427   if (distZ >= -halfTolerance && pzcut * vz >= << 562 
                                                   >> 563   C= B*B - 4.0*A*C;
                                                   >> 564   if (C > 0.0)
                                                   >> 565   {    
                                                   >> 566     G4double distR= (-B - std::sqrt(C)) / (2.0*A);
                                                   >> 567     G4double intZ = p.z()+distR*v.z();
                                                   >> 568     if ( (distR > halfRadTolerance)
                                                   >> 569       && (intZ >= zBottomCut-halfRadTolerance)
                                                   >> 570       && (intZ <= zTopCut+halfRadTolerance) )
                                                   >> 571     { 
                                                   >> 572       distMin = distR;
                                                   >> 573     }
                                                   >> 574     else if( (distR >- halfRadTolerance)
                                                   >> 575       && (intZ >= zBottomCut-halfRadTolerance)
                                                   >> 576       && (intZ <= zTopCut+halfRadTolerance) )
                                                   >> 577     {
                                                   >> 578       // p is on the curved surface, DistanceToIn returns 0 or kInfinity:
                                                   >> 579       // DistanceToIn returns 0, if second root is positive (means going inside)
                                                   >> 580       // If second root is negative, DistanceToIn returns kInfinity (outside)
                                                   >> 581       //
                                                   >> 582       distR = (-B + std::sqrt(C) ) / (2.0*A);
                                                   >> 583       if(distR>0.) { distMin=0.; }
                                                   >> 584     }
                                                   >> 585     else
                                                   >> 586     {
                                                   >> 587       distR= (-B + std::sqrt(C)) / (2.0*A);
                                                   >> 588       intZ = p.z()+distR*v.z();
                                                   >> 589       if ( (distR > halfRadTolerance)
                                                   >> 590         && (intZ >= zBottomCut-halfRadTolerance)
                                                   >> 591         && (intZ <= zTopCut+halfRadTolerance) )
                                                   >> 592       {
                                                   >> 593         G4ThreeVector norm=SurfaceNormal(p);
                                                   >> 594         if (norm.dot(v)<0.) { distMin = distR; }
                                                   >> 595       }
                                                   >> 596     }
                                                   >> 597     if ( (distMin!=kInfinity) && (distMin>dRmax) ) 
                                                   >> 598     {                    // Avoid rounding errors due to precision issues on
                                                   >> 599                          // 64 bits systems. Split long distances and recompute
                                                   >> 600       G4double fTerm = distMin-std::fmod(distMin,dRmax);
                                                   >> 601       distMin = fTerm + DistanceToIn(p+fTerm*v,v);
                                                   >> 602     }
                                                   >> 603   }
                                                   >> 604   
                                                   >> 605   if (std::fabs(distMin)<halfRadTolerance) { distMin=0.; }
                                                   >> 606   return distMin;
                                                   >> 607 } 
428                                                   608 
429   G4double rr = px * px + py * py + pz * pz;   << 609 ///////////////////////////////////////////////////////////////////////////////
430   G4double pv = px * vx + py * vy + pz * vz;   << 610 //
431   G4double distR = fQ1 * rr - fQ2;             << 611 // Calculate distance (<= actual) to closest surface of shape from outside
432   if (distR >= -halfTolerance && pv >= 0.) ret << 612 // - Return 0 if point inside
433                                                   613 
434   G4double A = vx * vx + vy * vy + vz * vz;    << 614 G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p) const
435   G4double B = pv;                             << 615 {
436   G4double C = rr - fR * fR;                   << 616   G4double distR, distZ;
437   G4double D = B * B - A * C;                  << 
438   // scratch^2 = R^2 - (R - halfTolerance)^2 = << 
439   G4double EPS = A * A * fR * kCarTolerance; / << 
440   if (D <= EPS) return kInfinity; // no inters << 
441                                                   617 
442   // Find intersection with Z planes           << 618   // normal vector:  parallel to normal, magnitude 1/(characteristic radius)
443   //                                              619   //
444   G4double invz  = (vz == 0) ? DBL_MAX : -1./v << 620   G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
445   G4double dz    = std::copysign(dzcut, invz); << 621                      p.y()/(ySemiAxis*ySemiAxis),
446   G4double tzmin = (pzcut - dz) * invz;        << 622                      p.z()/(zSemiAxis*zSemiAxis));
447   G4double tzmax = (pzcut + dz) * invz;        << 623   G4double radius= 1.0/norm.mag();
448                                                   624 
449   // Find intersection with lateral surface    << 625   // approximate distance to curved surface ( <= actual distance )
450   //                                              626   //
451   G4double tmp = -B - std::copysign(std::sqrt( << 627   distR= (p*norm - 1.0) * radius / 2.0;
452   G4double t1 = tmp / A;                       << 
453   G4double t2 = C / tmp;                       << 
454   G4double trmin = std::min(t1, t2);           << 
455   G4double trmax = std::max(t1, t2);           << 
456                                                   628 
457   // Return distance                           << 629   // Distance to z-cut plane
458   //                                              630   //
459   G4double tmin = std::max(tzmin, trmin);      << 631   distZ= zBottomCut - p.z();
460   G4double tmax = std::min(tzmax, trmax);      << 632   if (distZ < 0.0)
461                                                << 633   {
462   if (tmax - tmin <= halfTolerance) return kIn << 634     distZ = p.z() - zTopCut;
463   return (tmin < halfTolerance) ? offset : tmi << 635   }
464 }                                              << 
465                                                << 
466 ////////////////////////////////////////////// << 
467 //                                             << 
468 // Estimate distance to surface from outside   << 
469                                                   636 
470 G4double G4Ellipsoid::DistanceToIn(const G4Thr << 637   // Distance to closest surface from outside
471 {                                              << 638   //
472   G4double px = p.x();                         << 639   if (distZ < 0.0)
473   G4double py = p.y();                         << 640   {
474   G4double pz = p.z();                         << 641     return (distR < 0.0) ? 0.0 : distR;
475                                                << 642   }
476   // Safety distance to bounding box           << 643   else if (distR < 0.0)
477   G4double distX = std::abs(px) - fXmax;       << 644   {
478   G4double distY = std::abs(py) - fYmax;       << 645     return distZ;
479   G4double distZ = std::max(pz - fZTopCut, fZB << 646   }
480   G4double distB = std::max(std::max(distX, di << 647   else
481                                                << 648   {
482   // Safety distance to lateral surface        << 649     return (distZ < distR) ? distZ : distR;
483   G4double x = px * fSx;                       << 650   }
484   G4double y = py * fSy;                       << 
485   G4double z = pz * fSz;                       << 
486   G4double distR = std::sqrt(x*x + y*y + z*z)  << 
487                                                << 
488   // Return safety to in                       << 
489   G4double dist = std::max(distB, distR);      << 
490   return (dist < 0.) ? 0. : dist;              << 
491 }                                                 651 }
492                                                   652 
493 ////////////////////////////////////////////// << 653 ///////////////////////////////////////////////////////////////////////////////
494 //                                                654 //
495 // Calculate distance to surface from inside a << 655 // Calculate distance to surface of shape from `inside', allowing for tolerance
496                                                   656 
497 G4double G4Ellipsoid::DistanceToOut(const G4Th    657 G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p,
498                                     const G4Th    658                                     const G4ThreeVector& v,
499                                     const G4bo    659                                     const G4bool calcNorm,
500                                           G4bo << 660                                           G4bool *validNorm,
501                                           G4Th << 661                                           G4ThreeVector *n  ) const
502 {                                                 662 {
503   // Check if point flying away relative to Z  << 663   G4double distMin;
                                                   >> 664   enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
                                                   >> 665   
                                                   >> 666   distMin= kInfinity;
                                                   >> 667   surface= kNoSurf;
                                                   >> 668 
                                                   >> 669   // check to see if Z plane is relevant
504   //                                              670   //
505   G4double pz = p.z() * fSz;                   << 671   if (v.z() < 0.0)
506   G4double vz = v.z() * fSz;                   << 
507   G4double dzcut = fZDimCut;                   << 
508   G4double pzcut = pz - fZMidCut;              << 
509   G4double distZ = std::abs(pzcut) - dzcut;    << 
510   if (distZ >= -halfTolerance && pzcut * vz >  << 
511   {                                               672   {
512     if (calcNorm)                              << 673     G4double distZ = (zBottomCut - p.z()) / v.z();
                                                   >> 674     if (distZ < 0.0)
513     {                                             675     {
514       *validNorm = true;                       << 676       distZ= 0.0;
515       n->set(0., 0., std::copysign(1., pzcut)) << 677       if (!calcNorm) {return 0.0;}
516     }                                             678     }
517     return 0.;                                 << 679     distMin= distZ;
                                                   >> 680     surface= kPlaneSurf;
518   }                                               681   }
519                                                << 682   if (v.z() > 0.0)
520   // Check if point is flying away relative to << 
521   //                                           << 
522   G4double px = p.x() * fSx;                   << 
523   G4double py = p.y() * fSy;                   << 
524   G4double vx = v.x() * fSx;                   << 
525   G4double vy = v.y() * fSy;                   << 
526   G4double rr = px * px + py * py + pz * pz;   << 
527   G4double pv = px * vx + py * vy + pz * vz;   << 
528   G4double distR = fQ1 * rr - fQ2;             << 
529   if (distR >= -halfTolerance && pv > 0.)      << 
530   {                                               683   {
531     if (calcNorm)                              << 684     G4double distZ = (zTopCut - p.z()) / v.z();
                                                   >> 685     if (distZ < 0.0)
532     {                                             686     {
533       *validNorm = true;                       << 687       distZ= 0.0;
534       *n = G4ThreeVector(px*fSx, py*fSy, pz*fS << 688       if (!calcNorm) {return 0.0;}
535     }                                             689     }
536     return 0.;                                 << 690     distMin= distZ;
                                                   >> 691     surface= kPlaneSurf;
537   }                                               692   }
538                                                   693 
539   // Just in case check if point is outside (n << 694   // normal vector:  parallel to normal, magnitude 1/(characteristic radius)
540   //                                              695   //
541   if (std::max(distZ, distR) > halfTolerance)  << 696   G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis),
                                                   >> 697                          p.y()/(ySemiAxis*ySemiAxis),
                                                   >> 698                          p.z()/(zSemiAxis*zSemiAxis));
                                                   >> 699   
                                                   >> 700   // now check curved surface intercept
                                                   >> 701   //
                                                   >> 702   G4double A,B,C;
                                                   >> 703   
                                                   >> 704   A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
                                                   >> 705   C= (p * nearnorm) - 1.0;
                                                   >> 706   B= 2.0 * (v * nearnorm);
                                                   >> 707 
                                                   >> 708   C= B*B - 4.0*A*C;
                                                   >> 709   if (C > 0.0)
542   {                                               710   {
543 #ifdef G4SPECSDEBUG                            << 711     G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
544     std::ostringstream message;                << 712     if (distR < 0.0)
545     G4long oldprc = message.precision(16);     << 
546     message << "Point p is outside (!?) of sol << 
547             << GetName() << G4endl;            << 
548     message << "Position:  " << p << G4endl;;  << 
549     message << "Direction: " << v;             << 
550     G4cout.precision(oldprc);                  << 
551     G4Exception("G4Ellipsoid::DistanceToOut(p, << 
552                 JustWarning, message );        << 
553     DumpInfo();                                << 
554 #endif                                         << 
555     if (calcNorm)                              << 
556     {                                             713     {
557       *validNorm = true;                       << 714       distR= 0.0;
558       *n = ApproxSurfaceNormal(p);             << 715       if (!calcNorm) {return 0.0;}
                                                   >> 716     }
                                                   >> 717     if (distR < distMin)
                                                   >> 718     {
                                                   >> 719       distMin= distR;
                                                   >> 720       surface= kCurvedSurf;
559     }                                             721     }
560     return 0.;                                 << 
561   }                                               722   }
562                                                   723 
563   // Set coefficients of quadratic equation: A << 724   // set normal if requested
564   //                                              725   //
565   G4double A  = vx * vx + vy * vy + vz * vz;   << 726   if (calcNorm)
566   G4double B  = pv;                            << 
567   G4double C  = rr - fR * fR;                  << 
568   G4double D  = B * B - A * C;                 << 
569   // It is expected that the point is located  << 
570   // max term in the expression for discrimina << 
571   // max calculation error can be derived as f << 
572   // A * (1 + 2e) * R^2 * (1 + 2e) = A * R^2 + << 
573   G4double EPS = 4. * A * fR * fR * DBL_EPSILO << 
574                                                << 
575   if (D <= EPS) // no intersection             << 
576   {                                               727   {
577     if (calcNorm)                              << 728     if (surface == kNoSurf)
                                                   >> 729     {
                                                   >> 730       *validNorm = false;
                                                   >> 731     }
                                                   >> 732     else
578     {                                             733     {
579       *validNorm = true;                          734       *validNorm = true;
580       *n = G4ThreeVector(px*fSx, py*fSy, pz*fS << 735       switch (surface)
                                                   >> 736       {
                                                   >> 737         case kPlaneSurf:
                                                   >> 738           *n= G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
                                                   >> 739           break;
                                                   >> 740         case kCurvedSurf:
                                                   >> 741         {
                                                   >> 742           G4ThreeVector pexit= p + distMin*v;
                                                   >> 743           G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis),
                                                   >> 744                                  pexit.y()/(ySemiAxis*ySemiAxis),
                                                   >> 745                                  pexit.z()/(zSemiAxis*zSemiAxis));
                                                   >> 746           truenorm *= 1.0/truenorm.mag();
                                                   >> 747           *n= truenorm;
                                                   >> 748         } break;
                                                   >> 749         default:           // Should never reach this case ...
                                                   >> 750           DumpInfo();
                                                   >> 751           std::ostringstream message;
                                                   >> 752           G4int oldprc = message.precision(16);
                                                   >> 753           message << "Undefined side for valid surface normal to solid."
                                                   >> 754                   << G4endl
                                                   >> 755                   << "Position:"  << G4endl
                                                   >> 756                   << "   p.x() = "   << p.x()/mm << " mm" << G4endl
                                                   >> 757                   << "   p.y() = "   << p.y()/mm << " mm" << G4endl
                                                   >> 758                   << "   p.z() = "   << p.z()/mm << " mm" << G4endl
                                                   >> 759                   << "Direction:" << G4endl << G4endl
                                                   >> 760                   << "   v.x() = "   << v.x() << G4endl
                                                   >> 761                   << "   v.y() = "   << v.y() << G4endl
                                                   >> 762                   << "   v.z() = "   << v.z() << G4endl
                                                   >> 763                   << "Proposed distance :" << G4endl
                                                   >> 764                   << "   distMin = "    << distMin/mm << " mm";
                                                   >> 765           message.precision(oldprc);
                                                   >> 766           G4Exception("G4Ellipsoid::DistanceToOut(p,v,..)",
                                                   >> 767                       "GeomSolids1002", JustWarning, message);
                                                   >> 768           break;
                                                   >> 769       }
581     }                                             770     }
582     return 0.;                                 << 
583   }                                               771   }
                                                   >> 772    
                                                   >> 773   return distMin;
                                                   >> 774 }
                                                   >> 775 
                                                   >> 776 ///////////////////////////////////////////////////////////////////////////////
                                                   >> 777 //
                                                   >> 778 // Calculate distance (<=actual) to closest surface of shape from inside
                                                   >> 779 
                                                   >> 780 G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p) const
                                                   >> 781 {
                                                   >> 782   G4double distR, distZ;
                                                   >> 783 
                                                   >> 784 #ifdef G4SPECSDEBUG
                                                   >> 785   if( Inside(p) == kOutside )
                                                   >> 786   {
                                                   >> 787      DumpInfo();
                                                   >> 788      std::ostringstream message;
                                                   >> 789      G4int oldprc = message.precision(16);
                                                   >> 790      message << "Point p is outside !?" << G4endl
                                                   >> 791              << "Position:"  << G4endl
                                                   >> 792              << "   p.x() = "   << p.x()/mm << " mm" << G4endl
                                                   >> 793              << "   p.y() = "   << p.y()/mm << " mm" << G4endl
                                                   >> 794              << "   p.z() = "   << p.z()/mm << " mm";
                                                   >> 795      message.precision(oldprc) ;
                                                   >> 796      G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
                                                   >> 797                  JustWarning, message);
                                                   >> 798   }
                                                   >> 799 #endif
584                                                   800 
585   // Find intersection with Z cuts             << 801   // Normal vector:  parallel to normal, magnitude 1/(characteristic radius)
586   //                                              802   //
587   G4double tzmax = (vz == 0.) ? DBL_MAX : (std << 803   G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
                                                   >> 804                      p.y()/(ySemiAxis*ySemiAxis),
                                                   >> 805                      p.z()/(zSemiAxis*zSemiAxis));
588                                                   806 
589   // Find intersection with lateral surface    << 807   // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag())
590   //                                              808   //
591   G4double tmp = -B - std::copysign(std::sqrt( << 809   G4double radius= p.mag();
592   G4double trmax = (tmp < 0.) ? C/tmp : tmp/A; << 810   G4double tmp= norm.mag();
                                                   >> 811   if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
593                                                   812 
594   // Find distance and set normal, if required << 813   // Approximate distance to curved surface ( <= actual distance )
595   //                                              814   //
596   G4double tmax = std::min(tzmax, trmax);      << 815   distR = (1.0 - p*norm) * radius / 2.0;
597   //if (tmax < halfTolerance) tmax = 0.;       << 816     
                                                   >> 817   // Distance to z-cut plane
                                                   >> 818   //
                                                   >> 819   distZ = p.z() - zBottomCut;
                                                   >> 820   if (distZ < 0.0) {distZ= zTopCut - p.z();}
598                                                   821 
599   if (calcNorm)                                << 822   // Distance to closest surface from inside
                                                   >> 823   //
                                                   >> 824   if ( (distZ < 0.0) || (distR < 0.0) )
600   {                                               825   {
601     *validNorm = true;                         << 826     return 0.0;
602     if (tmax == tzmax)                         << 827   }
603     {                                          << 828   else
604       G4double pznew = pz + tmax * vz;         << 829   {
605       n->set(0., 0., (pznew > fZMidCut) ? 1. : << 830     return (distZ < distR) ? distZ : distR;
606     }                                          << 
607     else                                       << 
608     {                                          << 
609       G4double nx = (px + tmax * vx) * fSx;    << 
610       G4double ny = (py + tmax * vy) * fSy;    << 
611       G4double nz = (pz + tmax * vz) * fSz;    << 
612       *n = G4ThreeVector(nx, ny, nz).unit();   << 
613     }                                          << 
614   }                                               831   }
615   return tmax;                                 << 
616 }                                                 832 }
617                                                   833 
618 ////////////////////////////////////////////// << 834 ///////////////////////////////////////////////////////////////////////////////
619 //                                                835 //
620 // Estimate distance to surface from inside    << 836 // Create a List containing the transformed vertices
                                                   >> 837 // Ordering [0-3] -fDz cross section
                                                   >> 838 //          [4-7] +fDz cross section such that [0] is below [4],
                                                   >> 839 //                                             [1] below [5] etc.
                                                   >> 840 // Note:
                                                   >> 841 //  Caller has deletion resposibility
                                                   >> 842 //  Potential improvement: For last slice, use actual ending angle
                                                   >> 843 //                         to avoid rounding error problems.
                                                   >> 844 
                                                   >> 845 G4ThreeVectorList*
                                                   >> 846 G4Ellipsoid::CreateRotatedVertices(const G4AffineTransform& pTransform,
                                                   >> 847                                          G4int& noPolygonVertices) const
                                                   >> 848 {
                                                   >> 849   G4ThreeVectorList *vertices;
                                                   >> 850   G4ThreeVector vertex;
                                                   >> 851   G4double meshAnglePhi, meshRMaxFactor,
                                                   >> 852            crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
                                                   >> 853   G4double meshTheta, crossTheta, startTheta;
                                                   >> 854   G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
                                                   >> 855   G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
621                                                   856 
622 G4double G4Ellipsoid::DistanceToOut(const G4Th << 857   // Phi cross sections
623 {                                              << 858   //
624   // Safety distance in z direction            << 859   noPhiCrossSections=G4int (twopi/kMeshAngleDefault)+1;  // = 9!
625   G4double distZ = std::min(fZTopCut - p.z(),  << 860     
                                                   >> 861 /*
                                                   >> 862   if (noPhiCrossSections<kMinMeshSections)        // <3
                                                   >> 863   {
                                                   >> 864     noPhiCrossSections=kMinMeshSections;
                                                   >> 865   }
                                                   >> 866   else if (noPhiCrossSections>kMaxMeshSections)   // >37
                                                   >> 867   {
                                                   >> 868     noPhiCrossSections=kMaxMeshSections;
                                                   >> 869   }
                                                   >> 870 */
                                                   >> 871   meshAnglePhi=twopi/(noPhiCrossSections-1);
                                                   >> 872     
                                                   >> 873   // Set start angle such that mesh will be at fRMax
                                                   >> 874   // on the x axis. Will give better extent calculations when not rotated.
                                                   >> 875     
                                                   >> 876   sAnglePhi = -meshAnglePhi*0.5;
                                                   >> 877 
                                                   >> 878   // Theta cross sections
                                                   >> 879     
                                                   >> 880   noThetaSections = G4int(pi/kMeshAngleDefault)+3;  //  = 7!
                                                   >> 881 
                                                   >> 882 /*
                                                   >> 883   if (noThetaSections<kMinMeshSections)       // <3
                                                   >> 884   {
                                                   >> 885     noThetaSections=kMinMeshSections;
                                                   >> 886   }
                                                   >> 887   else if (noThetaSections>kMaxMeshSections)  // >37
                                                   >> 888   {
                                                   >> 889     noThetaSections=kMaxMeshSections;
                                                   >> 890   }
                                                   >> 891 */
                                                   >> 892   meshTheta= pi/(noThetaSections-2);
                                                   >> 893     
                                                   >> 894   // Set start angle such that mesh will be at fRMax
                                                   >> 895   // on the z axis. Will give better extent calculations when not rotated.
                                                   >> 896     
                                                   >> 897   startTheta = -meshTheta*0.5;
                                                   >> 898 
                                                   >> 899   meshRMaxFactor =  1.0/std::cos(0.5*
                                                   >> 900                     std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
                                                   >> 901   rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis);
                                                   >> 902   if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis;
                                                   >> 903   rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
                                                   >> 904   rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
                                                   >> 905   rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
                                                   >> 906   G4double* cosCrossTheta = new G4double[noThetaSections];
                                                   >> 907   G4double* sinCrossTheta = new G4double[noThetaSections];    
                                                   >> 908   vertices=new G4ThreeVectorList(noPhiCrossSections*noThetaSections);
                                                   >> 909   if (vertices && cosCrossTheta && sinCrossTheta)
                                                   >> 910   {
                                                   >> 911     for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
                                                   >> 912          crossSectionTheta++)
                                                   >> 913     {
                                                   >> 914       // Compute sine and cosine table (for historical reasons)
                                                   >> 915       //
                                                   >> 916       crossTheta=startTheta+crossSectionTheta*meshTheta;
                                                   >> 917       cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
                                                   >> 918       sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
                                                   >> 919     }
                                                   >> 920     for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
                                                   >> 921          crossSectionPhi++)
                                                   >> 922     {
                                                   >> 923       crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
                                                   >> 924       coscrossAnglePhi=std::cos(crossAnglePhi);
                                                   >> 925       sincrossAnglePhi=std::sin(crossAnglePhi);
                                                   >> 926       for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
                                                   >> 927            crossSectionTheta++)
                                                   >> 928       {
                                                   >> 929         // Compute coordinates of cross section at section crossSectionPhi
                                                   >> 930         //
                                                   >> 931         rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
                                                   >> 932         ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
                                                   >> 933         rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
                                                   >> 934         if (rz < zBottomCut)
                                                   >> 935           { rz= zBottomCut; }
                                                   >> 936         if (rz > zTopCut)
                                                   >> 937           { rz= zTopCut; }
                                                   >> 938         vertex= G4ThreeVector(rx,ry,rz);
                                                   >> 939         vertices->push_back(pTransform.TransformPoint(vertex));
                                                   >> 940       }    // Theta forward     
                                                   >> 941     }    // Phi
                                                   >> 942     noPolygonVertices = noThetaSections ;
                                                   >> 943   }
                                                   >> 944   else
                                                   >> 945   {
                                                   >> 946     DumpInfo();
                                                   >> 947     G4Exception("G4Ellipsoid::CreateRotatedVertices()",
                                                   >> 948                 "GeomSolids0003", FatalException,
                                                   >> 949                 "Error in allocation of vertices. Out of memory !");
                                                   >> 950   }
                                                   >> 951 
                                                   >> 952   delete[] cosCrossTheta;
                                                   >> 953   delete[] sinCrossTheta;
626                                                   954 
627   // Safety distance to lateral surface        << 955   return vertices;
628   G4double x = p.x() * fSx;                    << 
629   G4double y = p.y() * fSy;                    << 
630   G4double z = p.z() * fSz;                    << 
631   G4double distR = fR - std::sqrt(x*x + y*y +  << 
632                                                << 
633   // Return safety to out                      << 
634   G4double dist = std::min(distZ, distR);      << 
635   return (dist < 0.) ? 0. : dist;              << 
636 }                                                 956 }
637                                                   957 
638 //////////////////////////////////////////////    958 //////////////////////////////////////////////////////////////////////////
639 //                                                959 //
640 // Return entity type                          << 960 // G4EntityType
641                                                   961 
642 G4GeometryType G4Ellipsoid::GetEntityType() co    962 G4GeometryType G4Ellipsoid::GetEntityType() const
643 {                                                 963 {
644   return {"G4Ellipsoid"};                      << 964   return G4String("G4Ellipsoid");
645 }                                                 965 }
646                                                   966 
647 //////////////////////////////////////////////    967 //////////////////////////////////////////////////////////////////////////
648 //                                                968 //
649 // Make a clone of the object                     969 // Make a clone of the object
650                                                   970 
651 G4VSolid* G4Ellipsoid::Clone() const              971 G4VSolid* G4Ellipsoid::Clone() const
652 {                                                 972 {
653   return new G4Ellipsoid(*this);                  973   return new G4Ellipsoid(*this);
654 }                                                 974 }
655                                                   975 
656 //////////////////////////////////////////////    976 //////////////////////////////////////////////////////////////////////////
657 //                                                977 //
658 // Stream object contents to output stream     << 978 // Stream object contents to an output stream
659                                                   979 
660 std::ostream& G4Ellipsoid::StreamInfo( std::os    980 std::ostream& G4Ellipsoid::StreamInfo( std::ostream& os ) const
661 {                                                 981 {
662   G4long oldprc = os.precision(16);            << 982   G4int oldprc = os.precision(16);
663   os << "-------------------------------------    983   os << "-----------------------------------------------------------\n"
664      << "    *** Dump for solid - " << GetName    984      << "    *** Dump for solid - " << GetName() << " ***\n"
665      << "    =================================    985      << "    ===================================================\n"
666      << " Solid type: " << GetEntityType() <<  << 986      << " Solid type: G4Ellipsoid\n"
667      << " Parameters: \n"                         987      << " Parameters: \n"
668      << "    semi-axis x: " << GetDx()/mm << " << 988 
669      << "    semi-axis y: " << GetDy()/mm << " << 989      << "    semi-axis x: " << xSemiAxis/mm << " mm \n"
670      << "    semi-axis z: " << GetDz()/mm << " << 990      << "    semi-axis y: " << ySemiAxis/mm << " mm \n"
671      << "    lower cut in z: " << GetZBottomCu << 991      << "    semi-axis z: " << zSemiAxis/mm << " mm \n"
672      << "    upper cut in z: " << GetZTopCut() << 992      << "    max semi-axis: " << semiAxisMax/mm << " mm \n"
                                                   >> 993      << "    lower cut plane level z: " << zBottomCut/mm << " mm \n"
                                                   >> 994      << "    upper cut plane level z: " << zTopCut/mm << " mm \n"
673      << "-------------------------------------    995      << "-----------------------------------------------------------\n";
674   os.precision(oldprc);                           996   os.precision(oldprc);
                                                   >> 997 
675   return os;                                      998   return os;
676 }                                                 999 }
677                                                   1000 
678 ////////////////////////////////////////////// << 1001 ////////////////////////////////////////////////////////////////////
679 //                                                1002 //
680 // Return volume                               << 1003 // GetPointOnSurface
681                                                   1004 
682 G4double G4Ellipsoid::GetCubicVolume()         << 1005 G4ThreeVector G4Ellipsoid::GetPointOnSurface() const
683 {                                                 1006 {
684   if (fCubicVolume == 0.)                      << 1007   G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
                                                   >> 1008   G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3;
                                                   >> 1009 
                                                   >> 1010   max1  = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
                                                   >> 1011   max1  = max1 > zSemiAxis ? max1 : zSemiAxis;
                                                   >> 1012   if (max1 == xSemiAxis)      { max2 = ySemiAxis; max3 = zSemiAxis; }
                                                   >> 1013   else if (max1 == ySemiAxis) { max2 = xSemiAxis; max3 = zSemiAxis; }
                                                   >> 1014   else                        { max2 = xSemiAxis; max3 = ySemiAxis; }
                                                   >> 1015 
                                                   >> 1016   phi   = G4RandFlat::shoot(0.,twopi);
                                                   >> 1017   
                                                   >> 1018   cosphi = std::cos(phi);   sinphi = std::sin(phi);
                                                   >> 1019   costheta = G4RandFlat::shoot(zBottomCut,zTopCut)/zSemiAxis;
                                                   >> 1020   sintheta = std::sqrt(1.-sqr(costheta));
                                                   >> 1021   
                                                   >> 1022   alpha = 1.-sqr(max2/max1); beta  = 1.-sqr(max3/max1);
                                                   >> 1023   
                                                   >> 1024   aTop    = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis));
                                                   >> 1025   aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis));
                                                   >> 1026   
                                                   >> 1027   // approximation
                                                   >> 1028   // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf"
                                                   >> 1029   aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)-
                                                   >> 1030                             1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta)));
                                                   >> 1031 
                                                   >> 1032   aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis);
                                                   >> 1033   
                                                   >> 1034   if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis )
                                                   >> 1035    || ( zTopCut == 0 && zBottomCut ==0 ) )
685   {                                               1036   {
686     G4double piAB_3 = CLHEP::pi * fDx * fDy /  << 1037     aTop = 0; aBottom = 0;
687     fCubicVolume = 4. * piAB_3 * fDz;          << 
688     if (fZBottomCut > -fDz)                    << 
689     {                                          << 
690       G4double hbot = 1. + fZBottomCut / fDz;  << 
691       fCubicVolume -= piAB_3 * hbot * hbot * ( << 
692     }                                          << 
693     if (fZTopCut < fDz)                        << 
694     {                                          << 
695       G4double htop = 1. - fZTopCut / fDz;     << 
696       fCubicVolume -= piAB_3 * htop * htop * ( << 
697     }                                          << 
698   }                                               1038   }
699   return fCubicVolume;                         << 1039   
700 }                                              << 1040   chose = G4RandFlat::shoot(0.,aTop + aBottom + aCurved); 
701                                                << 1041   
702 ////////////////////////////////////////////// << 1042   if(chose < aCurved)
703 //                                             << 1043   { 
704 // Calculate area of lateral surface           << 1044     xRand = xSemiAxis*sintheta*cosphi;
705                                                << 1045     yRand = ySemiAxis*sintheta*sinphi;
706 G4double G4Ellipsoid::LateralSurfaceArea() con << 1046     zRand = zSemiAxis*costheta;
707 {                                              << 1047     return G4ThreeVector (xRand,yRand,zRand); 
708   constexpr G4int NPHI = 1000.;                << 
709   constexpr G4double dPhi = CLHEP::halfpi/NPHI << 
710   constexpr G4double eps = 4.*DBL_EPSILON;     << 
711                                                << 
712   G4double aa = fDx*fDx;                       << 
713   G4double bb = fDy*fDy;                       << 
714   G4double cc = fDz*fDz;                       << 
715   G4double ab = fDx*fDy;                       << 
716   G4double cc_aa = cc/aa;                      << 
717   G4double cc_bb = cc/bb;                      << 
718   G4double zmax = std::min(fZTopCut, fDz);     << 
719   G4double zmin = std::max(fZBottomCut,-fDz);  << 
720   G4double zmax_c = zmax/fDz;                  << 
721   G4double zmin_c = zmin/fDz;                  << 
722   G4double area = 0.;                          << 
723                                                << 
724   if (aa == bb) // spheroid, use analytical ex << 
725   {                                            << 
726     G4double k = fDz/fDx;                      << 
727     G4double kk = k*k;                         << 
728     if (kk < 1. - eps)                         << 
729     {                                          << 
730       G4double invk = fDx/fDz;                 << 
731       G4double root = std::sqrt(1. - kk);      << 
732       G4double tmax = zmax_c*root;             << 
733       G4double tmin = zmin_c*root;             << 
734       area = CLHEP::pi*ab*                     << 
735         ((zmax_c*std::sqrt(kk + tmax*tmax) - z << 
736          (std::asinh(tmax*invk) - std::asinh(t << 
737     }                                          << 
738     else if (kk > 1. + eps)                    << 
739     {                                          << 
740       G4double invk = fDx/fDz;                 << 
741       G4double root = std::sqrt(kk - 1.);      << 
742       G4double tmax = zmax_c*root;             << 
743       G4double tmin = zmin_c*root;             << 
744       area = CLHEP::pi*ab*                     << 
745         ((zmax_c*std::sqrt(kk - tmax*tmax) - z << 
746          (std::asin(tmax*invk) - std::asin(tmi << 
747     }                                          << 
748     else                                       << 
749     {                                          << 
750       area = CLHEP::twopi*fDx*(zmax - zmin);   << 
751     }                                          << 
752     return area;                               << 
753   }                                               1048   }
754                                                << 1049   else if(chose >= aCurved && chose < aCurved + aTop)
755   // ellipsoid, integration along phi          << 
756   for (G4int i = 0; i < NPHI; ++i)             << 
757   {                                               1050   {
758     G4double sinPhi = std::sin(dPhi*(i + 0.5)) << 1051     xRand = G4RandFlat::shoot(-1.,1.)*xSemiAxis
759     G4double kk = cc_aa + (cc_bb - cc_aa)*sinP << 1052           * std::sqrt(1-sqr(zTopCut/zSemiAxis));
760     if (kk < 1. - eps)                         << 1053     yRand = G4RandFlat::shoot(-1.,1.)*ySemiAxis
761     {                                          << 1054           * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis));
762       G4double root = std::sqrt(1. - kk);      << 1055     zRand = zTopCut;
763       G4double tmax = zmax_c*root;             << 1056     return G4ThreeVector (xRand,yRand,zRand);
764       G4double tmin = zmin_c*root;             << 
765       G4double invk = 1./std::sqrt(kk);        << 
766       area += 2.*ab*dPhi*                      << 
767         ((zmax_c*std::sqrt(kk + tmax*tmax) - z << 
768          (std::asinh(tmax*invk) - std::asinh(t << 
769     }                                          << 
770     else if (kk > 1. + eps)                    << 
771     {                                          << 
772       G4double root = std::sqrt(kk - 1.);      << 
773       G4double tmax = zmax_c*root;             << 
774       G4double tmin = zmin_c*root;             << 
775       G4double invk = 1./std::sqrt(kk);        << 
776       area += 2.*ab*dPhi*                      << 
777         ((zmax_c*std::sqrt(kk - tmax*tmax) - z << 
778          (std::asin(tmax*invk) - std::asin(tmi << 
779     }                                          << 
780     else                                       << 
781     {                                          << 
782       area += 4.*ab*dPhi*(zmax_c - zmin_c);    << 
783     }                                          << 
784   }                                               1057   }
785   return area;                                 << 1058   else
786 }                                              << 
787                                                << 
788 ////////////////////////////////////////////// << 
789 //                                             << 
790 // Return surface area                         << 
791                                                << 
792 G4double G4Ellipsoid::GetSurfaceArea()         << 
793 {                                              << 
794   if (fSurfaceArea == 0.)                      << 
795   {                                               1059   {
796     G4double piAB = CLHEP::pi * fDx * fDy;     << 1060     xRand = G4RandFlat::shoot(-1.,1.)*xSemiAxis
797     fSurfaceArea = LateralSurfaceArea();       << 1061           * std::sqrt(1-sqr(zBottomCut/zSemiAxis));
798     if (fZBottomCut > -fDz)                    << 1062     yRand = G4RandFlat::shoot(-1.,1.)*ySemiAxis
799     {                                          << 1063           * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis)); 
800       G4double hbot = 1. + fZBottomCut / fDz;  << 1064     zRand = zBottomCut;
801       fSurfaceArea += piAB * hbot * (2. - hbot << 1065     return G4ThreeVector (xRand,yRand,zRand);
802     }                                          << 
803     if (fZTopCut < fDz)                        << 
804     {                                          << 
805       G4double htop = 1. - fZTopCut / fDz;     << 
806       fSurfaceArea += piAB * htop * (2. - htop << 
807     }                                          << 
808   }                                               1066   }
809   return fSurfaceArea;                         << 
810 }                                                 1067 }
811                                                   1068 
812 ////////////////////////////////////////////// << 1069 /////////////////////////////////////////////////////////////////////////////
813 //                                             << 
814 // Return random point on surface              << 
815                                                << 
816 G4ThreeVector G4Ellipsoid::GetPointOnSurface() << 
817 {                                              << 
818   G4double A    = GetDx();                     << 
819   G4double B    = GetDy();                     << 
820   G4double C    = GetDz();                     << 
821   G4double Zbot = GetZBottomCut();             << 
822   G4double Ztop = GetZTopCut();                << 
823                                                << 
824   // Calculate cut areas                       << 
825   G4double Hbot = 1. + Zbot / C;               << 
826   G4double Htop = 1. - Ztop / C;               << 
827   G4double piAB = CLHEP::pi * A * B;           << 
828   G4double Sbot = piAB * Hbot * (2. - Hbot);   << 
829   G4double Stop = piAB * Htop * (2. - Htop);   << 
830                                                << 
831   // Get area of lateral surface               << 
832   if (fLateralArea == 0.)                      << 
833   {                                            << 
834     G4AutoLock l(&lateralareaMutex);           << 
835     fLateralArea = LateralSurfaceArea();       << 
836     l.unlock();                                << 
837   }                                            << 
838   G4double Slat = fLateralArea;                << 
839                                                << 
840   // Select surface (0 - bottom cut, 1 - later << 
841   G4double select = (Sbot + Slat + Stop) * G4Q << 
842   G4int k = 0;                                 << 
843   if (select > Sbot) k = 1;                    << 
844   if (select > Sbot + Slat) k = 2;             << 
845                                                << 
846   // Pick random point on selected surface (re << 
847   G4ThreeVector p;                             << 
848   switch (k)                                   << 
849   {                                            << 
850     case 0: // bootom z-cut                    << 
851     {                                          << 
852       G4double scale = std::sqrt(Hbot * (2. -  << 
853       G4TwoVector rho = G4RandomPointInEllipse << 
854       p.set(rho.x(), rho.y(), Zbot);           << 
855       break;                                   << 
856     }                                          << 
857     case 1: // lateral surface                 << 
858     {                                          << 
859       G4double x, y, z;                        << 
860       G4double mu_max = std::max(std::max(A *  << 
861       for (G4int i = 0; i < 1000; ++i)         << 
862       {                                        << 
863         // generate random point on unit spher << 
864         z = (Zbot + (Ztop - Zbot) * G4QuickRan << 
865         G4double rho = std::sqrt((1. + z) * (1 << 
866         G4double phi = CLHEP::twopi * G4QuickR << 
867         x = rho * std::cos(phi);               << 
868         y = rho * std::sin(phi);               << 
869         // check  acceptance                   << 
870         G4double xbc = x * B * C;              << 
871         G4double yac = y * A * C;              << 
872         G4double zab = z * A * B;              << 
873         G4double mu  = std::sqrt(xbc * xbc + y << 
874         if (mu_max * G4QuickRand() <= mu) brea << 
875       }                                        << 
876       p.set(A * x, B * y, C * z);              << 
877       break;                                   << 
878     }                                          << 
879     case 2: // top z-cut                       << 
880     {                                          << 
881       G4double scale  = std::sqrt(Htop * (2. - << 
882       G4TwoVector rho = G4RandomPointInEllipse << 
883       p.set(rho.x(), rho.y(), Ztop);           << 
884       break;                                   << 
885     }                                          << 
886   }                                            << 
887   return p;                                    << 
888 }                                              << 
889                                                << 
890 ////////////////////////////////////////////// << 
891 //                                                1070 //
892 // Methods for visualisation                      1071 // Methods for visualisation
893                                                   1072 
894 void G4Ellipsoid::DescribeYourselfTo (G4VGraph    1073 void G4Ellipsoid::DescribeYourselfTo (G4VGraphicsScene& scene) const
895 {                                                 1074 {
896   scene.AddSolid(*this);                          1075   scene.AddSolid(*this);
897 }                                                 1076 }
898                                                   1077 
899 ////////////////////////////////////////////// << 
900 //                                             << 
901 // Return vis extent                           << 
902                                                << 
903 G4VisExtent G4Ellipsoid::GetExtent() const        1078 G4VisExtent G4Ellipsoid::GetExtent() const
904 {                                                 1079 {
905   return { -fXmax, fXmax, -fYmax, fYmax, fZBot << 1080   // Define the sides of the box into which the G4Ellipsoid instance would fit.
                                                   >> 1081   //
                                                   >> 1082   return G4VisExtent (-semiAxisMax, semiAxisMax,
                                                   >> 1083                       -semiAxisMax, semiAxisMax,
                                                   >> 1084                       -semiAxisMax, semiAxisMax);
906 }                                                 1085 }
907                                                   1086 
908 ////////////////////////////////////////////// << 
909 //                                             << 
910 // Create polyhedron                           << 
911                                                << 
912 G4Polyhedron* G4Ellipsoid::CreatePolyhedron ()    1087 G4Polyhedron* G4Ellipsoid::CreatePolyhedron () const
913 {                                                 1088 {
914   return new G4PolyhedronEllipsoid(fDx, fDy, f << 1089   return new G4PolyhedronEllipsoid(xSemiAxis, ySemiAxis, zSemiAxis,
                                                   >> 1090                                    zBottomCut, zTopCut);
915 }                                                 1091 }
916                                                   1092 
917 ////////////////////////////////////////////// << 
918 //                                             << 
919 // Return pointer to polyhedron                << 
920                                                << 
921 G4Polyhedron* G4Ellipsoid::GetPolyhedron () co    1093 G4Polyhedron* G4Ellipsoid::GetPolyhedron () const
922 {                                                 1094 {
923   if (fpPolyhedron == nullptr ||               << 1095   if (!fpPolyhedron ||
924       fRebuildPolyhedron ||                       1096       fRebuildPolyhedron ||
925       fpPolyhedron->GetNumberOfRotationStepsAt    1097       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
926       fpPolyhedron->GetNumberOfRotationSteps()    1098       fpPolyhedron->GetNumberOfRotationSteps())
927     {                                             1099     {
928       G4AutoLock l(&polyhedronMutex);             1100       G4AutoLock l(&polyhedronMutex);
929       delete fpPolyhedron;                        1101       delete fpPolyhedron;
930       fpPolyhedron = CreatePolyhedron();          1102       fpPolyhedron = CreatePolyhedron();
931       fRebuildPolyhedron = false;                 1103       fRebuildPolyhedron = false;
932       l.unlock();                                 1104       l.unlock();
933     }                                             1105     }
934   return fpPolyhedron;                            1106   return fpPolyhedron;
935 }                                                 1107 }
936                                                << 
937 #endif // !defined(G4GEOM_USE_UELLIPSOID) || ! << 
938                                                   1108