Geant4 Cross Reference

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


  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 // G4EllipticalTube implementation             << 
 27 //                                                 26 //
 28 // Author: David C. Williams (davidw@scipp.ucs <<  27 // $Id: G4EllipticalTube.cc,v 1.26 2006/06/29 18:48:31 gunter Exp $
 29 // Revision: Evgueni Tcherniaev (evgueni.tcher <<  28 // GEANT4 tag $Name: geant4-08-01 $
                                                   >>  29 //
                                                   >>  30 // 
                                                   >>  31 // --------------------------------------------------------------------
                                                   >>  32 // GEANT 4 class source file
                                                   >>  33 //
                                                   >>  34 //
                                                   >>  35 // G4EllipticalTube.cc
                                                   >>  36 //
                                                   >>  37 // Implementation of a CSG volume representing a tube with elliptical cross
                                                   >>  38 // section (geant3 solid 'ELTU')
                                                   >>  39 //
 30 // -------------------------------------------     40 // --------------------------------------------------------------------
 31                                                    41 
 32 #include "G4EllipticalTube.hh"                     42 #include "G4EllipticalTube.hh"
 33                                                    43 
 34 #if !(defined(G4GEOM_USE_UELLIPTICALTUBE) && d << 
 35                                                << 
 36 #include "G4GeomTools.hh"                      << 
 37 #include "G4RandomTools.hh"                    << 
 38 #include "G4ClippablePolygon.hh"                   44 #include "G4ClippablePolygon.hh"
 39 #include "G4AffineTransform.hh"                    45 #include "G4AffineTransform.hh"
                                                   >>  46 #include "G4SolidExtentList.hh"
 40 #include "G4VoxelLimits.hh"                        47 #include "G4VoxelLimits.hh"
 41 #include "G4BoundingEnvelope.hh"               <<  48 #include "meshdefs.hh"
 42                                                    49 
 43 #include "Randomize.hh"                            50 #include "Randomize.hh"
 44                                                    51 
 45 #include "G4VGraphicsScene.hh"                     52 #include "G4VGraphicsScene.hh"
                                                   >>  53 #include "G4Polyhedron.hh"
 46 #include "G4VisExtent.hh"                          54 #include "G4VisExtent.hh"
 47                                                    55 
 48 #include "G4AutoLock.hh"                       << 
 49                                                << 
 50 namespace                                      << 
 51 {                                              << 
 52   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE << 
 53 }                                              << 
 54                                                << 
 55 using namespace CLHEP;                             56 using namespace CLHEP;
 56                                                    57 
 57 ////////////////////////////////////////////// << 
 58 //                                                 58 //
 59 // Constructor                                     59 // Constructor
 60                                                <<  60 //
 61 G4EllipticalTube::G4EllipticalTube( const G4St <<  61 G4EllipticalTube::G4EllipticalTube( const G4String &name, 
 62                                           G4do <<  62                                           G4double theDx,
 63                                           G4do <<  63                                           G4double theDy,
 64                                           G4do <<  64                                           G4double theDz )
 65   : G4VSolid(name), fDx(Dx), fDy(Dy), fDz(Dz)  <<  65   : G4VSolid( name ), fCubicVolume(0.), fpPolyhedron(0)
 66 {                                              <<  66 {
 67   CheckParameters();                           <<  67   dx = theDx;
                                                   >>  68   dy = theDy;
                                                   >>  69   dz = theDz;
 68 }                                                  70 }
 69                                                    71 
 70 ////////////////////////////////////////////// <<  72 
 71 //                                                 73 //
 72 // Fake default constructor - sets only member     74 // Fake default constructor - sets only member data and allocates memory
 73 //                            for usage restri     75 //                            for usage restricted to object persistency.
 74                                                <<  76 //
 75 G4EllipticalTube::G4EllipticalTube( __void__&      77 G4EllipticalTube::G4EllipticalTube( __void__& a )
 76   : G4VSolid(a), halfTolerance(0.), fDx(0.), f <<  78   : G4VSolid(a), fCubicVolume(0.), fpPolyhedron(0)
 77     fRsph(0.), fDDx(0.), fDDy(0.), fSx(0.), fS << 
 78     fQ1(0.), fQ2(0.), fScratch(0.)             << 
 79 {                                                  79 {
 80 }                                                  80 }
 81                                                    81 
 82 ////////////////////////////////////////////// << 
 83 //                                                 82 //
 84 // Destructor                                      83 // Destructor
 85                                                << 
 86 G4EllipticalTube::~G4EllipticalTube()          << 
 87 {                                              << 
 88   delete fpPolyhedron; fpPolyhedron = nullptr; << 
 89 }                                              << 
 90                                                << 
 91 ////////////////////////////////////////////// << 
 92 //                                                 84 //
 93 // Copy constructor                            <<  85 G4EllipticalTube::~G4EllipticalTube() {delete fpPolyhedron;}
 94                                                    86 
 95 G4EllipticalTube::G4EllipticalTube(const G4Ell << 
 96   : G4VSolid(rhs), halfTolerance(rhs.halfToler << 
 97     fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),  << 
 98     fCubicVolume(rhs.fCubicVolume), fSurfaceAr << 
 99     fRsph(rhs.fRsph), fDDx(rhs.fDDx), fDDy(rhs << 
100     fSx(rhs.fSx), fSy(rhs.fSy), fR(rhs.fR),    << 
101     fQ1(rhs.fQ1), fQ2(rhs.fQ2), fScratch(rhs.f << 
102 {                                              << 
103 }                                              << 
104                                                    87 
105 ////////////////////////////////////////////// << 
106 //                                                 88 //
107 // Assignment operator                         <<  89 // CalculateExtent
108                                                << 
109 G4EllipticalTube& G4EllipticalTube::operator = << 
110 {                                              << 
111    // Check assignment to self                 << 
112    //                                          << 
113    if (this == &rhs)  { return *this; }        << 
114                                                << 
115    // Copy base class data                     << 
116    //                                          << 
117    G4VSolid::operator=(rhs);                   << 
118                                                << 
119    // Copy data                                << 
120    //                                          << 
121    halfTolerance = rhs.halfTolerance;          << 
122    fDx = rhs.fDx;                              << 
123    fDy = rhs.fDy;                              << 
124    fDz = rhs.fDz;                              << 
125    fCubicVolume = rhs.fCubicVolume;            << 
126    fSurfaceArea = rhs.fSurfaceArea;            << 
127                                                << 
128    fRsph = rhs.fRsph;                          << 
129    fDDx  = rhs.fDDx;                           << 
130    fDDy  = rhs.fDDy;                           << 
131    fSx   = rhs.fSx;                            << 
132    fSy   = rhs.fSy;                            << 
133    fR    = rhs.fR;                             << 
134    fQ1   = rhs.fQ1;                            << 
135    fQ2   = rhs.fQ2;                            << 
136    fScratch = rhs.fScratch;                    << 
137                                                << 
138    fRebuildPolyhedron = false;                 << 
139    delete fpPolyhedron; fpPolyhedron = nullptr << 
140                                                << 
141    return *this;                               << 
142 }                                              << 
143                                                << 
144 ////////////////////////////////////////////// << 
145 //                                                 90 //
146 // Check dimensions                            <<  91 G4bool
                                                   >>  92 G4EllipticalTube::CalculateExtent( const EAxis axis,
                                                   >>  93                                    const G4VoxelLimits &voxelLimit,
                                                   >>  94                                    const G4AffineTransform &transform,
                                                   >>  95                                          G4double &min, G4double &max ) const
                                                   >>  96 {
                                                   >>  97   G4SolidExtentList  extentList( axis, voxelLimit );
                                                   >>  98   
                                                   >>  99   //
                                                   >> 100   // We are going to divide up our elliptical face into small
                                                   >> 101   // pieces
                                                   >> 102   //
                                                   >> 103   
                                                   >> 104   //
                                                   >> 105   // Choose phi size of our segment(s) based on constants as
                                                   >> 106   // defined in meshdefs.hh
                                                   >> 107   //
                                                   >> 108   G4int numPhi = kMaxMeshSections;
                                                   >> 109   G4double sigPhi = twopi/numPhi;
                                                   >> 110   
                                                   >> 111   //
                                                   >> 112   // We have to be careful to keep our segments completely outside
                                                   >> 113   // of the elliptical surface. To do so we imagine we have
                                                   >> 114   // a simple (unit radius) circular cross section (as in G4Tubs) 
                                                   >> 115   // and then "stretch" the dimensions as necessary to fit the ellipse.
                                                   >> 116   //
                                                   >> 117   G4double rFudge = 1.0/std::cos(0.5*sigPhi);
                                                   >> 118   G4double dxFudge = dx*rFudge,
                                                   >> 119            dyFudge = dy*rFudge;
                                                   >> 120   
                                                   >> 121   //
                                                   >> 122   // As we work around the elliptical surface, we build
                                                   >> 123   // a "phi" segment on the way, and keep track of two
                                                   >> 124   // additional polygons for the two ends.
                                                   >> 125   //
                                                   >> 126   G4ClippablePolygon endPoly1, endPoly2, phiPoly;
                                                   >> 127   
                                                   >> 128   G4double phi = 0, 
                                                   >> 129            cosPhi = std::cos(phi),
                                                   >> 130            sinPhi = std::sin(phi);
                                                   >> 131   G4ThreeVector v0( dxFudge*cosPhi, dyFudge*sinPhi, +dz ),
                                                   >> 132                 v1( dxFudge*cosPhi, dyFudge*sinPhi, -dz ),
                                                   >> 133                 w0, w1;
                                                   >> 134   transform.ApplyPointTransform( v0 );
                                                   >> 135   transform.ApplyPointTransform( v1 );
                                                   >> 136   do
                                                   >> 137   {
                                                   >> 138     phi += sigPhi;
                                                   >> 139     if (numPhi == 1) phi = 0;  // Try to avoid roundoff
                                                   >> 140     cosPhi = std::cos(phi), 
                                                   >> 141     sinPhi = std::sin(phi);
                                                   >> 142     
                                                   >> 143     w0 = G4ThreeVector( dxFudge*cosPhi, dyFudge*sinPhi, +dz );
                                                   >> 144     w1 = G4ThreeVector( dxFudge*cosPhi, dyFudge*sinPhi, -dz );
                                                   >> 145     transform.ApplyPointTransform( w0 );
                                                   >> 146     transform.ApplyPointTransform( w1 );
                                                   >> 147     
                                                   >> 148     //
                                                   >> 149     // Add a point to our z ends
                                                   >> 150     //
                                                   >> 151     endPoly1.AddVertexInOrder( v0 );
                                                   >> 152     endPoly2.AddVertexInOrder( v1 );
                                                   >> 153     
                                                   >> 154     //
                                                   >> 155     // Build phi polygon
                                                   >> 156     //
                                                   >> 157     phiPoly.ClearAllVertices();
                                                   >> 158     
                                                   >> 159     phiPoly.AddVertexInOrder( v0 );
                                                   >> 160     phiPoly.AddVertexInOrder( v1 );
                                                   >> 161     phiPoly.AddVertexInOrder( w1 );
                                                   >> 162     phiPoly.AddVertexInOrder( w0 );
                                                   >> 163     
                                                   >> 164     if (phiPoly.PartialClip( voxelLimit, axis ))
                                                   >> 165     {
                                                   >> 166       //
                                                   >> 167       // Get unit normal
                                                   >> 168       //
                                                   >> 169       phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
                                                   >> 170       
                                                   >> 171       extentList.AddSurface( phiPoly );
                                                   >> 172     }
                                                   >> 173 
                                                   >> 174     //
                                                   >> 175     // Next vertex
                                                   >> 176     //    
                                                   >> 177     v0 = w0;
                                                   >> 178     v1 = w1;
                                                   >> 179   } while( --numPhi > 0 );
147                                                   180 
148 void G4EllipticalTube::CheckParameters()       << 
149 {                                              << 
150   // Check dimensions                          << 
151   //                                              181   //
152   halfTolerance = 0.5*kCarTolerance; // half t << 182   // Process the end pieces
153   G4double dmin = 2*kCarTolerance;             << 183   //
154   if (fDx < dmin || fDy < dmin || fDz < dmin)  << 184   if (endPoly1.PartialClip( voxelLimit, axis ))
155   {                                               185   {
156     std::ostringstream message;                << 186     static const G4ThreeVector normal(0,0,+1);
157     message << "Invalid (too small or negative << 187     endPoly1.SetNormal( transform.TransformAxis(normal) );
158             << GetName()                       << 188     extentList.AddSurface( endPoly1 );
159             << "\n  Dx = " << fDx              << 
160             << "\n  Dy = " << fDy              << 
161             << "\n  Dz = " << fDz;             << 
162     G4Exception("G4EllipticalTube::CheckParame << 
163           FatalException, message);            << 
164   }                                               189   }
165                                                << 190   
166   // Set pre-calculatated values               << 191   if (endPoly2.PartialClip( voxelLimit, axis ))
                                                   >> 192   {
                                                   >> 193     static const G4ThreeVector normal(0,0,-1);
                                                   >> 194     endPoly2.SetNormal( transform.TransformAxis(normal) );
                                                   >> 195     extentList.AddSurface( endPoly2 );
                                                   >> 196   }
                                                   >> 197   
167   //                                              198   //
168   halfTolerance = 0.5*kCarTolerance; // half t << 199   // Return min/max value
169   fRsph = std::sqrt(fDx * fDx + fDy * fDy + fD << 200   //
170   fDDx = fDx * fDx; // X semi-axis squared     << 201   return extentList.GetExtent( min, max );
171   fDDy = fDy * fDy; // Y semi-axis squared     << 
172                                                << 
173   fR = std::min(fDx, fDy); // resulting radius << 
174   fSx = fR / fDx; // X scale factor            << 
175   fSy = fR / fDy; // Y scale factor            << 
176                                                << 
177   fQ1 = 0.5 / fR; // distance approxiamtion di << 
178   fQ2 = 0.5 * (fR + halfTolerance * halfTolera << 
179   fScratch = 2. * fR * fR * DBL_EPSILON; // sc << 
180   // fScratch = (B * B / A) * (2. + halfTolera << 
181 }                                                 202 }
182                                                   203 
183 ////////////////////////////////////////////// << 
184 //                                             << 
185 // Get bounding box                            << 
186                                                   204 
187 void G4EllipticalTube::BoundingLimits( G4Three << 205 //
188                                        G4Three << 206 // Inside
                                                   >> 207 //
                                                   >> 208 // Note that for this solid, we've decided to define the tolerant
                                                   >> 209 // surface as that which is bounded by ellipses with axes
                                                   >> 210 // at +/- 0.5*kCarTolerance.
                                                   >> 211 //
                                                   >> 212 EInside G4EllipticalTube::Inside( const G4ThreeVector& p ) const
189 {                                                 213 {
190   pMin.set(-fDx,-fDy,-fDz);                    << 214   static const G4double halfTol = 0.5*kCarTolerance;
191   pMax.set( fDx, fDy, fDz);                    << 215   
                                                   >> 216   //
                                                   >> 217   // Check z extents: are we outside?
                                                   >> 218   //
                                                   >> 219   G4double absZ = std::fabs(p.z());
                                                   >> 220   if (absZ > dz+halfTol) return kOutside;
                                                   >> 221   
                                                   >> 222   //
                                                   >> 223   // Check x,y: are we outside?
                                                   >> 224   //
                                                   >> 225   // G4double x = p.x(), y = p.y();
                                                   >> 226   
                                                   >> 227   if (CheckXY(p.x(), p.y(), +halfTol) > 1.0) return kOutside;
                                                   >> 228   
                                                   >> 229   //
                                                   >> 230   // We are either inside or on the surface: recheck z extents
                                                   >> 231   //
                                                   >> 232   if (absZ > dz-halfTol) return kSurface;
                                                   >> 233   
                                                   >> 234   //
                                                   >> 235   // Recheck x,y
                                                   >> 236   //
                                                   >> 237   if (CheckXY(p.x(), p.y(), -halfTol) > 1.0) return kSurface;
                                                   >> 238   
                                                   >> 239   return kInside;
192 }                                                 240 }
193                                                   241 
194 ////////////////////////////////////////////// << 
195 //                                             << 
196 // Calculate extent under transform and specif << 
197                                                   242 
198 G4bool                                         << 243 //
199 G4EllipticalTube::CalculateExtent( const EAxis << 244 // SurfaceNormal
200                                    const G4Vox << 245 //
201                                    const G4Aff << 246 G4ThreeVector G4EllipticalTube::SurfaceNormal( const G4ThreeVector& p ) const
202                                          G4dou << 
203 {                                                 247 {
204   G4ThreeVector bmin, bmax;                    << 
205   G4bool exist;                                << 
206                                                << 
207   // Check bounding box (bbox)                 << 
208   //                                              248   //
209   BoundingLimits(bmin,bmax);                   << 249   // Which of the three surfaces are we closest to (approximately)?
210   G4BoundingEnvelope bbox(bmin,bmax);          << 
211 #ifdef G4BBOX_EXTENT                           << 
212   return bbox.CalculateExtent(pAxis,pVoxelLimi << 
213 #endif                                         << 
214   if (bbox.BoundingBoxVsVoxelLimits(pAxis, pVo << 
215   {                                            << 
216     return exist = pMin < pMax;                << 
217   }                                            << 
218                                                << 
219   G4double dx = fDx;                           << 
220   G4double dy = fDy;                           << 
221   G4double dz = fDz;                           << 
222                                                << 
223   // Set bounding envelope (benv) and calculat << 
224   //                                              250   //
225   const G4int NSTEPS = 24; // number of steps  << 251   G4double distZ = std::fabs(p.z()) - dz;
226   G4double ang = twopi/NSTEPS;                 << 252   
227                                                << 253   G4double rxy = CheckXY( p.x(), p.y() );
228   G4double sinHalf = std::sin(0.5*ang);        << 254   G4double distR2 = (rxy < DBL_MIN) ? DBL_MAX : 1.0/rxy;
229   G4double cosHalf = std::cos(0.5*ang);        << 
230   G4double sinStep = 2.*sinHalf*cosHalf;       << 
231   G4double cosStep = 1. - 2.*sinHalf*sinHalf;  << 
232   G4double sx = dx/cosHalf;                    << 
233   G4double sy = dy/cosHalf;                    << 
234                                                << 
235   G4double sinCur = sinHalf;                   << 
236   G4double cosCur = cosHalf;                   << 
237   G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS << 
238   for (G4int k=0; k<NSTEPS; ++k)               << 
239   {                                            << 
240     baseA[k].set(sx*cosCur,sy*sinCur,-dz);     << 
241     baseB[k].set(sx*cosCur,sy*sinCur, dz);     << 
242                                                   255 
243     G4double sinTmp = sinCur;                  << 256   //
244     sinCur = sinCur*cosStep + cosCur*sinStep;  << 257   // Closer to z?
245     cosCur = cosCur*cosStep - sinTmp*sinStep;  << 258   //
246   }                                            << 259   if (distZ*distZ < distR2) 
                                                   >> 260     return G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
247                                                   261 
248   std::vector<const G4ThreeVectorList *> polyg << 262   //
249   polygons[0] = &baseA;                        << 263   // Closer to x/y
250   polygons[1] = &baseB;                        << 264   //
251   G4BoundingEnvelope benv(bmin, bmax, polygons << 265   return G4ThreeVector( p.x()*dy*dy, p.y()*dx*dx, 0.0 ).unit();
252   exist = benv.CalculateExtent(pAxis, pVoxelLi << 
253   return exist;                                << 
254 }                                                 266 }
255                                                   267 
256 ////////////////////////////////////////////// << 268 
257 //                                                269 //
258 // Determine where is point: inside, outside o << 270 // DistanceToIn(p,v)
259 //                                                271 //
260                                                << 272 // Unlike DistanceToOut(p,v), it is possible for the trajectory
261 EInside G4EllipticalTube::Inside( const G4Thre << 273 // to miss. The geometric calculations here are quite simple.
262 {                                              << 274 // More difficult is the logic required to prevent particles
263   G4double x = p.x() * fSx;                    << 275 // from sneaking (or leaking) between the elliptical and end
264   G4double y = p.y() * fSy;                    << 276 // surfaces.
265   G4double distR = fQ1 * (x * x + y * y) - fQ2 << 
266   G4double distZ = std::abs(p.z()) - fDz;      << 
267   G4double dist = std::max(distR, distZ);      << 
268                                                << 
269   if (dist > halfTolerance) return kOutside;   << 
270   return (dist > -halfTolerance) ? kSurface :  << 
271 }                                              << 
272                                                << 
273 ////////////////////////////////////////////// << 
274 //                                                277 //
275 // Return unit normal at surface closest to p  << 278 // Keep in mind that the true distance is allowed to be
276                                                << 279 // negative if the point is currently on the surface. For oblique
277 G4ThreeVector G4EllipticalTube::SurfaceNormal( << 280 // angles, it can be very negative. 
                                                   >> 281 //
                                                   >> 282 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p,
                                                   >> 283                                          const G4ThreeVector& v ) const
278 {                                                 284 {
279   G4ThreeVector norm(0, 0, 0);                 << 285   static const G4double halfTol = 0.5*kCarTolerance;
280   G4int nsurf = 0;                             << 286     
                                                   >> 287   //
                                                   >> 288   // Check z = -dz planer surface
                                                   >> 289   //
                                                   >> 290   G4double sigz = p.z()+dz;
281                                                   291 
282   // check lateral surface                     << 292   if (sigz < halfTol)
283   G4double x = p.x() * fSx;                    << 
284   G4double y = p.y() * fSy;                    << 
285   G4double distR = fQ1 * (x * x + y * y) - fQ2 << 
286   if (std::abs(distR) <= halfTolerance)        << 
287   {                                               293   {
288     norm = G4ThreeVector(p.x() * fDDy, p.y() * << 294     //
289     ++nsurf;                                   << 295     // We are "behind" the shape in z, and so can
                                                   >> 296     // potentially hit the rear face. Correct direction?
                                                   >> 297     //
                                                   >> 298     if (v.z() <= 0)
                                                   >> 299     {
                                                   >> 300       //
                                                   >> 301       // As long as we are far enough away, we know we
                                                   >> 302       // can't intersect
                                                   >> 303       //
                                                   >> 304       if (sigz < 0) return kInfinity;
                                                   >> 305       
                                                   >> 306       //
                                                   >> 307       // Otherwise, we don't intersect unless we are
                                                   >> 308       // on the surface of the ellipse
                                                   >> 309       //
                                                   >> 310       if (CheckXY(p.x(),p.y(),-halfTol) <= 1.0) return kInfinity;
                                                   >> 311     }
                                                   >> 312     else
                                                   >> 313     {
                                                   >> 314       //
                                                   >> 315       // How far?
                                                   >> 316       //
                                                   >> 317       G4double s = -sigz/v.z();
                                                   >> 318       
                                                   >> 319       //
                                                   >> 320       // Where does that place us?
                                                   >> 321       //
                                                   >> 322       G4double xi = p.x() + s*v.x(),
                                                   >> 323                yi = p.y() + s*v.y();
                                                   >> 324       
                                                   >> 325       //
                                                   >> 326       // Is this on the surface (within ellipse)?
                                                   >> 327       //
                                                   >> 328       if (CheckXY(xi,yi) <= 1.0)
                                                   >> 329       {
                                                   >> 330         //
                                                   >> 331         // Yup. Return s, unless we are on the surface
                                                   >> 332         //
                                                   >> 333         return (sigz < -halfTol) ? s : 0;
                                                   >> 334       }
                                                   >> 335       else if (xi*dy*dy*v.x() + yi*dx*dx*v.y() >= 0)
                                                   >> 336       {
                                                   >> 337         //
                                                   >> 338         // Else, if we are traveling outwards, we know
                                                   >> 339         // we must miss
                                                   >> 340         //
                                                   >> 341         return kInfinity;
                                                   >> 342       }
                                                   >> 343     }
290   }                                               344   }
291                                                   345 
292   // check lateral bases                       << 346   //
293   G4double distZ = std::abs(p.z()) - fDz;      << 347   // Check z = +dz planer surface
294   if (std::abs(distZ) <= halfTolerance)        << 348   //
                                                   >> 349   sigz = p.z() - dz;
                                                   >> 350   
                                                   >> 351   if (sigz > -halfTol)
295   {                                               352   {
296     norm.setZ(p.z() < 0 ? -1. : 1.);           << 353     if (v.z() >= 0)
297     ++nsurf;                                   << 354     {
                                                   >> 355       if (sigz > 0) return kInfinity;
                                                   >> 356       if (CheckXY(p.x(),p.y(),-halfTol) <= 1.0) return kInfinity;
                                                   >> 357     }
                                                   >> 358     else {
                                                   >> 359       G4double s = -sigz/v.z();
                                                   >> 360 
                                                   >> 361       G4double xi = p.x() + s*v.x(),
                                                   >> 362                yi = p.y() + s*v.y();
                                                   >> 363       
                                                   >> 364       if (CheckXY(xi,yi) <= 1.0)
                                                   >> 365       {
                                                   >> 366         return (sigz > -halfTol) ? s : 0;
                                                   >> 367       }
                                                   >> 368       else if (xi*dy*dy*v.x() + yi*dx*dx*v.y() >= 0)
                                                   >> 369       {
                                                   >> 370         return kInfinity;
                                                   >> 371       }
                                                   >> 372     }
298   }                                               373   }
                                                   >> 374   
                                                   >> 375   //
                                                   >> 376   // Check intersection with the elliptical tube
                                                   >> 377   //
                                                   >> 378   G4double s[2];
                                                   >> 379   G4int n = IntersectXY( p, v, s );
                                                   >> 380   
                                                   >> 381   if (n==0) return kInfinity;
                                                   >> 382   
                                                   >> 383   //
                                                   >> 384   // Is the original point on the surface?
                                                   >> 385   //
                                                   >> 386   if (std::fabs(p.z()) < dz+halfTol) {
                                                   >> 387     if (CheckXY( p.x(), p.y(), halfTol ) < 1.0)
                                                   >> 388     {
                                                   >> 389       //
                                                   >> 390       // Well, yes, but are we traveling inwards at this point?
                                                   >> 391       //
                                                   >> 392       if (p.x()*dy*dy*v.x() + p.y()*dx*dx*v.y() < 0) return 0;
                                                   >> 393     }
                                                   >> 394   }
                                                   >> 395   
                                                   >> 396   //
                                                   >> 397   // We are now certain that point p is not on the surface of 
                                                   >> 398   // the solid (and thus std::fabs(s[0]) > halfTol). 
                                                   >> 399   // Return kInfinity if the intersection is "behind" the point.
                                                   >> 400   //
                                                   >> 401   if (s[0] < 0) return kInfinity;
                                                   >> 402   
                                                   >> 403   //
                                                   >> 404   // Check to see if we intersect the tube within
                                                   >> 405   // dz, but only when we know it might miss
                                                   >> 406   //
                                                   >> 407   G4double zi = p.z() + s[0]*v.z();
299                                                   408 
300   // return normal                             << 409   if (v.z() < 0)
301   if (nsurf == 1) return norm;                 << 
302   else if (nsurf > 1) return norm.unit(); // e << 
303   else                                         << 
304   {                                               410   {
305     // Point is not on the surface             << 411     if (zi < -dz) return kInfinity;
306     //                                         << 412   }
307 #ifdef G4SPECDEBUG                             << 413   else if (v.z() > 0)
308     std::ostringstream message;                << 414   {
309     G4long oldprc = message.precision(16);     << 415     if (zi > +dz) return kInfinity;
310     message << "Point p is not on surface (!?) << 
311             << GetName() << G4endl;            << 
312     message << "Position:\n";                  << 
313     message << "   p.x() = " << p.x()/mm << "  << 
314     message << "   p.y() = " << p.y()/mm << "  << 
315     message << "   p.z() = " << p.z()/mm << "  << 
316     G4cout.precision(oldprc);                  << 
317     G4Exception("G4EllipticalTube::SurfaceNorm << 
318                 JustWarning, message );        << 
319     DumpInfo();                                << 
320 #endif                                         << 
321     return ApproxSurfaceNormal(p);             << 
322   }                                               416   }
323 }                                              << 
324                                                   417 
325 ////////////////////////////////////////////// << 418   return s[0];
326 //                                             << 
327 // Find surface nearest to point and return co << 
328 // The algorithm is similar to the algorithm u << 
329 // This method normally should not be called.  << 
330                                                << 
331 G4ThreeVector                                  << 
332 G4EllipticalTube::ApproxSurfaceNormal( const G << 
333 {                                              << 
334   G4double x = p.x() * fSx;                    << 
335   G4double y = p.y() * fSy;                    << 
336   G4double distR = fQ1 * (x * x + y * y) - fQ2 << 
337   G4double distZ = std::abs(p.z()) - fDz;      << 
338   if (distR > distZ && (x * x + y * y) > 0)    << 
339     return G4ThreeVector(p.x() * fDDy, p.y() * << 
340   else                                         << 
341     return {0, 0, (p.z() < 0 ? -1. : 1.)};     << 
342 }                                                 419 }
343                                                   420 
344 ////////////////////////////////////////////// << 
345 //                                             << 
346 // Calculate distance to shape from outside, a << 
347 // return kInfinity if no intersection, or dis << 
348                                                   421 
349 G4double G4EllipticalTube::DistanceToIn( const << 422 //
350                                          const << 423 // DistanceToIn(p)
                                                   >> 424 //
                                                   >> 425 // The distance from a point to an ellipse (in 2 dimensions) is a
                                                   >> 426 // surprisingly complicated quadric expression (this is easy to
                                                   >> 427 // appreciate once one understands that there may be up to
                                                   >> 428 // four lines normal to the ellipse intersecting any point). To 
                                                   >> 429 // solve it exactly would be rather time consuming. This method, 
                                                   >> 430 // however, is supposed to be a quick check, and is allowed to be an
                                                   >> 431 // underestimate.
                                                   >> 432 //
                                                   >> 433 // So, I will use the following underestimate of the distance
                                                   >> 434 // from an outside point to an ellipse. First: find the intersection "A"
                                                   >> 435 // of the line from the origin to the point with the ellipse.
                                                   >> 436 // Find the line passing through "A" and tangent to the ellipse 
                                                   >> 437 // at A. The distance of the point p from the ellipse will be approximated
                                                   >> 438 // as the distance to this line.
                                                   >> 439 //
                                                   >> 440 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p ) const
351 {                                                 441 {
352   G4double offset = 0.;                        << 442   static const G4double halfTol = 0.5*kCarTolerance;
353   G4ThreeVector pcur = p;                      << 443   
354                                                << 444   if (CheckXY( p.x(), p.y(), +halfTol ) < 1.0)
355   // Check if point is flying away             << 
356   //                                           << 
357   G4double safex = std::abs(pcur.x()) - fDx;   << 
358   G4double safey = std::abs(pcur.y()) - fDy;   << 
359   G4double safez = std::abs(pcur.z()) - fDz;   << 
360                                                << 
361   if (safez >= -halfTolerance && pcur.z() * v. << 
362   if (safey >= -halfTolerance && pcur.y() * v. << 
363   if (safex >= -halfTolerance && pcur.x() * v. << 
364                                                << 
365   // Relocate point, if required               << 
366   //                                           << 
367   G4double Dmax = 32. * fRsph;                 << 
368   if (std::max(std::max(safex, safey), safez)  << 
369   {                                               445   {
370     offset = (1. - 1.e-08) * pcur.mag() - 2. * << 446     //
371     pcur += offset * v;                        << 447     // We are inside or on the surface of the
372     G4double dist = DistanceToIn(pcur, v);     << 448     // elliptical cross section in x/y. Check z
373     return (dist == kInfinity) ? kInfinity : d << 449     //
                                                   >> 450     if (p.z() < -dz-halfTol) 
                                                   >> 451       return -p.z()-dz;
                                                   >> 452     else if (p.z() > dz+halfTol)
                                                   >> 453       return p.z()-dz;
                                                   >> 454     else
                                                   >> 455       return 0;    // On any surface here (or inside)
374   }                                               456   }
375                                                << 457   
376   // Scale elliptical tube to cylinder         << 
377   //                                              458   //
378   G4double px = pcur.x() * fSx;                << 459   // Find point on ellipse
379   G4double py = pcur.y() * fSy;                << 
380   G4double pz = pcur.z();                      << 
381   G4double vx = v.x() * fSx;                   << 
382   G4double vy = v.y() * fSy;                   << 
383   G4double vz = v.z();                         << 
384                                                << 
385   // Set coefficients of quadratic equation: A << 
386   //                                              460   //
387   G4double rr = px * px + py * py;             << 461   G4double qnorm = CheckXY( p.x(), p.y() );
388   G4double A  = vx * vx + vy * vy;             << 462   if (qnorm < DBL_MIN) return 0;  // This should never happen
389   G4double B  = px * vx + py * vy;             << 463   
390   G4double C  = rr - fR * fR;                  << 464   G4double q = 1.0/std::sqrt(qnorm);
391   G4double D  = B * B - A * C;                 << 465   
392                                                << 466   G4double xe = q*p.x(), ye = q*p.y();
393   // Check if point is flying away relative to << 467      
394   //                                              468   //
395   G4double distR  = fQ1 * rr - fQ2;            << 469   // Get tangent to ellipse
396   G4bool parallelToZ = (A < DBL_EPSILON || std << 
397   if (distR >= -halfTolerance && (B >= 0. || p << 
398                                                << 
399   // Find intersection with Z planes           << 
400   //                                              470   //
401   G4double invz  = (vz == 0) ? DBL_MAX : -1./v << 471   G4double tx = -ye*dx*dx, ty = +xe*dy*dy;
402   G4double dz    = std::copysign(fDz, invz);   << 472   G4double tnorm = std::sqrt( tx*tx + ty*ty );
403   G4double tzmin = (pz - dz) * invz;           << 473   
404   G4double tzmax = (pz + dz) * invz;           << 
405                                                << 
406   // Solve qudratic equation. There are two ca << 
407   //   1) trajectory parallel to Z axis (A = 0 << 
408   //   2) touch (D = 0) or no intersection (D  << 
409   //                                              474   //
410   if (parallelToZ) return (tzmin<halfTolerance << 475   // Calculate distance
411   if (D <= A * A * fScratch) return kInfinity; << 476   //
412                                                << 477   G4double distR = ( (p.x()-xe)*ty - (p.y()-ye)*tx )/tnorm;
413   // Find roots of quadratic equation          << 478   
414   G4double tmp = -B - std::copysign(std::sqrt( << 479   //
415   G4double t1 = tmp / A;                       << 480   // Add the result in quadrature if we are, in addition,
416   G4double t2 = C / tmp;                       << 481   // outside the z bounds of the shape
417   G4double trmin = std::min(t1, t2);           << 482   //
418   G4double trmax = std::max(t1, t2);           << 483   // We could save some time by returning the maximum rather
419                                                << 484   // than the quadrature sum
420   // Return distance                           << 485   //
421   G4double tin  = std::max(tzmin, trmin);      << 486   if (p.z() < -dz) 
422   G4double tout = std::min(tzmax, trmax);      << 487     return std::sqrt( (p.z()+dz)*(p.z()+dz) + distR*distR );
                                                   >> 488   else if (p.z() > dz)
                                                   >> 489     return std::sqrt( (p.z()-dz)*(p.z()-dz) + distR*distR );
423                                                   490 
424   if (tout <= tin + halfTolerance) return kInf << 491   return distR;
425   return (tin<halfTolerance) ? offset : tin +  << 
426 }                                                 492 }
427                                                   493 
428 ////////////////////////////////////////////// << 
429 //                                             << 
430 // Estimate distance to the surface from outsi << 
431 // returns 0 if point is inside                << 
432                                                << 
433 G4double G4EllipticalTube::DistanceToIn( const << 
434 {                                              << 
435   // safety distance to bounding box           << 
436   G4double distX = std::abs(p.x()) - fDx;      << 
437   G4double distY = std::abs(p.y()) - fDy;      << 
438   G4double distZ = std::abs(p.z()) - fDz;      << 
439   G4double distB = std::max(std::max(distX, di << 
440   // return (distB < 0) ? 0 : distB;           << 
441                                                << 
442   // safety distance to lateral surface        << 
443   G4double x = p.x() * fSx;                    << 
444   G4double y = p.y() * fSy;                    << 
445   G4double distR = std::sqrt(x * x + y * y) -  << 
446                                                << 
447   // return SafetyToIn                         << 
448   G4double dist = std::max(distB, distR);      << 
449   return (dist < 0) ? 0 : dist;                << 
450 }                                              << 
451                                                << 
452 ////////////////////////////////////////////// << 
453 //                                             << 
454 // Calculate distance to shape from inside and << 
455 // at exit point, if required                  << 
456 // - when leaving the surface, return 0        << 
457                                                   494 
                                                   >> 495 //
                                                   >> 496 // DistanceToOut(p,v)
                                                   >> 497 //
                                                   >> 498 // This method can be somewhat complicated for a general shape.
                                                   >> 499 // For a convex one, like this, there are several simplifications,
                                                   >> 500 // the most important of which is that one can treat the surfaces
                                                   >> 501 // as infinite in extent when deciding if the p is on the surface.
                                                   >> 502 //
458 G4double G4EllipticalTube::DistanceToOut( cons    503 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p,
459                                           cons    504                                           const G4ThreeVector& v,
460                                           cons    505                                           const G4bool calcNorm,
461                                                << 506                                                 G4bool *validNorm,
462                                                << 507                                                 G4ThreeVector *norm ) const
463 {                                                 508 {
464   // Check if point flying away relative to Z  << 509   static const G4double halfTol = 0.5*kCarTolerance;
                                                   >> 510   
465   //                                              511   //
466   G4double pz = p.z();                         << 512   // Our normal is always valid
467   G4double vz = v.z();                         << 
468   G4double distZ = std::abs(pz) - fDz;         << 
469   if (distZ >= -halfTolerance && pz * vz > 0)  << 
470   {                                            << 
471     if (calcNorm)                              << 
472     {                                          << 
473       *validNorm = true;                       << 
474       n->set(0, 0, (pz < 0) ? -1. : 1.);       << 
475     }                                          << 
476     return 0.;                                 << 
477   }                                            << 
478   G4double tzmax = (vz == 0) ? DBL_MAX : (std: << 
479                                                << 
480   // Scale elliptical tube to cylinder         << 
481   //                                              513   //
482   G4double px = p.x() * fSx;                   << 514   if (calcNorm) *validNorm = true;
483   G4double py = p.y() * fSy;                   << 515   
484   G4double vx = v.x() * fSx;                   << 516   G4double sBest = kInfinity;
485   G4double vy = v.y() * fSy;                   << 517   const G4ThreeVector *nBest=0;
486                                                << 518   
487   // Check if point is flying away relative to << 
488   //                                              519   //
489   G4double rr = px * px + py * py;             << 520   // Might we intersect the -dz surface?
490   G4double B  = px * vx + py * vy;             << 521   //
491   G4double distR  = fQ1 * rr - fQ2;            << 522   if (v.z() < 0)
492   if (distR >= -halfTolerance && B > 0.)       << 
493   {                                               523   {
494     if (calcNorm)                              << 524     static const G4ThreeVector normHere(0.0,0.0,-1.0);
495     {                                          << 525     //
496       *validNorm = true;                       << 526     // Yup. What distance?
497       *n = G4ThreeVector(px * fDDy, py * fDDx, << 527     //
                                                   >> 528     sBest = -(p.z()+dz)/v.z();
                                                   >> 529     
                                                   >> 530     //
                                                   >> 531     // Are we on the surface? If so, return zero
                                                   >> 532     //
                                                   >> 533     if (p.z() < -dz+halfTol) {
                                                   >> 534       if (calcNorm) *norm = normHere;
                                                   >> 535       return 0;
498     }                                             536     }
499     return 0.;                                 << 537     else
                                                   >> 538       nBest = &normHere;
500   }                                               539   }
501                                                << 540   
502   // Just in case check if point is outside, n << 
503   //                                              541   //
504   if (std::max(distZ, distR) > halfTolerance)  << 542   // How about the +dz surface?
                                                   >> 543   //
                                                   >> 544   if (v.z() > 0)
505   {                                               545   {
506 #ifdef G4SPECDEBUG                             << 546     static const G4ThreeVector normHere(0.0,0.0,+1.0);
507     std::ostringstream message;                << 547     //
508     G4long oldprc = message.precision(16);     << 548     // Yup. What distance?
509     message << "Point p is outside (!?) of sol << 549     //
510             << GetName() << G4endl;            << 550     G4double s = (dz-p.z())/v.z();
511     message << "Position:  " << p << G4endl;;  << 551     
512     message << "Direction: " << v;             << 552     //
513     G4cout.precision(oldprc);                  << 553     // Are we on the surface? If so, return zero
514     G4Exception("G4EllipticalTube::DistanceToO << 554     //
515                 JustWarning, message );        << 555     if (p.z() > +dz-halfTol) {
516     DumpInfo();                                << 556       if (calcNorm) *norm = normHere;
517 #endif                                         << 557       return 0;
518     if (calcNorm)                              << 
519     {                                          << 
520       *validNorm = true;                       << 
521       *n = ApproxSurfaceNormal(p);             << 
522     }                                             558     }
523     return 0.;                                 << 559     
                                                   >> 560     //
                                                   >> 561     // Best so far?
                                                   >> 562     //
                                                   >> 563     if (s < sBest) { sBest = s; nBest = &normHere; }
524   }                                               564   }
525                                                << 565   
526   // Set coefficients of quadratic equation: A << 
527   //                                              566   //
528   G4double A  = vx * vx + vy * vy;             << 567   // Check furthest intersection with ellipse 
529   G4double C  = rr - fR * fR;                  << 
530   G4double D  = B * B - A * C;                 << 
531                                                << 
532   // Solve qudratic equation. There are two sp << 
533   //   1) trajectory parallel to Z axis (A = 0 << 
534   //   2) touch (D = 0) or no intersection (D  << 
535   //                                              568   //
536   G4bool parallelToZ = (A < DBL_EPSILON || std << 569   G4double s[2];
537   if (parallelToZ) // 1)                       << 570   G4int n = IntersectXY( p, v, s );
                                                   >> 571 
                                                   >> 572   if (n == 0)
538   {                                               573   {
539     if (calcNorm)                              << 574     if (sBest == kInfinity)
540     {                                             575     {
541       *validNorm = true;                       << 576       G4cout.precision(16) ;
542       n->set(0, 0, (vz < 0) ? -1. : 1.);       << 577       G4cout << G4endl ;
                                                   >> 578       DumpInfo();
                                                   >> 579       G4cout << "Position:"  << G4endl << G4endl ;
                                                   >> 580       G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
                                                   >> 581       G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
                                                   >> 582       G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
                                                   >> 583       G4cout << "Direction:" << G4endl << G4endl;
                                                   >> 584       G4cout << "v.x() = "   << v.x() << G4endl;
                                                   >> 585       G4cout << "v.y() = "   << v.y() << G4endl;
                                                   >> 586       G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
                                                   >> 587       G4cout << "Proposed distance :" << G4endl << G4endl;
                                                   >> 588       G4cout << "snxt = "    << sBest/mm << " mm" << G4endl << G4endl;
                                                   >> 589       G4Exception( "G4EllipticalTube::DistanceToOut(p,v,...)",
                                                   >> 590                    "Notification", JustWarning, "Point p is outside !?" );
543     }                                             591     }
544     return tzmax;                              << 592     if (calcNorm) *norm = *nBest;
                                                   >> 593     return sBest;
545   }                                               594   }
546   if (D <= A * A * fScratch) // 2)             << 595   else if (s[n-1] > sBest)
547   {                                               596   {
548     if (calcNorm)                              << 597     if (calcNorm) *norm = *nBest;
549     {                                          << 598     return sBest;
550       *validNorm = true;                       << 599   }  
551       *n = G4ThreeVector(px * fDDy, py * fDDx, << 600   sBest = s[n-1];
552     }                                          << 601       
553     return 0.;                                 << 602   //
554   }                                            << 603   // Intersection with ellipse. Get normal at intersection point.
555                                                << 
556   // Find roots of quadratic equation          << 
557   G4double tmp = -B - std::copysign(std::sqrt( << 
558   G4double t1 = tmp / A;                       << 
559   G4double t2 = C / tmp;                       << 
560   G4double trmax = std::max(t1, t2);           << 
561                                                << 
562   // Return distance                           << 
563   G4double tmax = std::min(tzmax, trmax);      << 
564                                                << 
565   // Set normal, if required, and return dista << 
566   //                                              604   //
567   if (calcNorm)                                   605   if (calcNorm)
568   {                                               606   {
569     *validNorm = true;                         << 607     G4ThreeVector ip = p + sBest*v;
570     G4ThreeVector pnew = p + tmax * v;         << 608     *norm = G4ThreeVector( ip.x()*dy*dy, ip.y()*dx*dx, 0.0 ).unit();
571     if (tmax == tzmax)                         << 609   }
572       n->set(0, 0, (pnew.z() < 0) ? -1. : 1.); << 610   
573     else                                       << 611   //
574       *n = G4ThreeVector(pnew.x() * fDDy, pnew << 612   // Do we start on the surface?
                                                   >> 613   //
                                                   >> 614   if (CheckXY( p.x(), p.y(), -halfTol ) > 1.0)
                                                   >> 615   {
                                                   >> 616     //
                                                   >> 617     // Well, yes, but are we traveling outwards at this point?
                                                   >> 618     //
                                                   >> 619     if (p.x()*dy*dy*v.x() + p.y()*dx*dx*v.y() > 0) return 0;
575   }                                               620   }
576   return tmax;                                 << 621   
                                                   >> 622   return sBest;
577 }                                                 623 }
578                                                   624 
579 ////////////////////////////////////////////// << 625 
580 //                                                626 //
581 // Estimate distance to the surface from insid << 627 // DistanceToOut(p)
582 // returns 0 if point is outside               << 628 //
                                                   >> 629 // See DistanceToIn(p) for notes on the distance from a point
                                                   >> 630 // to an ellipse in two dimensions.
                                                   >> 631 //
                                                   >> 632 // The approximation used here for a point inside the ellipse
                                                   >> 633 // is to find the intersection with the ellipse of the lines 
                                                   >> 634 // through the point and parallel to the x and y axes. The
                                                   >> 635 // distance of the point from the line connecting the two 
                                                   >> 636 // intersecting points is then used.
583 //                                                637 //
584                                                << 
585 G4double G4EllipticalTube::DistanceToOut( cons    638 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p ) const
586 {                                                 639 {
587 #ifdef G4SPECDEBUG                             << 640   static const G4double halfTol = 0.5*kCarTolerance;
588   if( Inside(p) == kOutside )                  << 641   
                                                   >> 642   //
                                                   >> 643   // We need to calculate the distances to all surfaces,
                                                   >> 644   // and then return the smallest
                                                   >> 645   //
                                                   >> 646   // Check -dz and +dz surface
                                                   >> 647   //
                                                   >> 648   G4double sBest = dz - std::fabs(p.z());
                                                   >> 649   if (sBest < halfTol) return 0;
                                                   >> 650   
                                                   >> 651   //
                                                   >> 652   // Check elliptical surface: find intersection of
                                                   >> 653   // line through p and parallel to x axis
                                                   >> 654   //
                                                   >> 655   G4double radical = 1.0 - p.y()*p.y()/dy/dy;
                                                   >> 656   if (radical < +DBL_MIN) return 0;
                                                   >> 657   
                                                   >> 658   G4double xi = dx*std::sqrt( radical );
                                                   >> 659   if (p.x() < 0) xi = -xi;
                                                   >> 660   
                                                   >> 661   //
                                                   >> 662   // Do the same with y axis
                                                   >> 663   //
                                                   >> 664   radical = 1.0 - p.x()*p.x()/dx/dx;
                                                   >> 665   if (radical < +DBL_MIN) return 0;
                                                   >> 666   
                                                   >> 667   G4double yi = dy*std::sqrt( radical );
                                                   >> 668   if (p.y() < 0) yi = -yi;
                                                   >> 669   
                                                   >> 670   //
                                                   >> 671   // Get distance from p to the line connecting
                                                   >> 672   // these two points
                                                   >> 673   //
                                                   >> 674   G4double xdi = p.x() - xi,
                                                   >> 675      ydi = yi - p.y();
                                                   >> 676 
                                                   >> 677   G4double normi = std::sqrt( xdi*xdi + ydi*ydi );
                                                   >> 678   if (normi < halfTol) return 0;
                                                   >> 679   xdi /= normi;
                                                   >> 680   ydi /= normi;
                                                   >> 681   
                                                   >> 682   G4double s = 0.5*(xdi*(p.y()-yi) - ydi*(p.x()-xi));
                                                   >> 683   if (xi*yi < 0) s = -s;
                                                   >> 684   
                                                   >> 685   if (s < sBest) sBest = s;
                                                   >> 686   
                                                   >> 687   //
                                                   >> 688   // Return best answer
                                                   >> 689   //
                                                   >> 690   return sBest < halfTol ? 0 : sBest;
                                                   >> 691 }
                                                   >> 692 
                                                   >> 693 
                                                   >> 694 //
                                                   >> 695 // IntersectXY
                                                   >> 696 //
                                                   >> 697 // Decide if and where the x/y trajectory hits the elliptical cross
                                                   >> 698 // section.
                                                   >> 699 //
                                                   >> 700 // Arguments:
                                                   >> 701 //     p     - (in) Point on trajectory
                                                   >> 702 //     v     - (in) Vector along trajectory
                                                   >> 703 //     s     - (out) Up to two points of intersection, where the
                                                   >> 704 //                   intersection point is p + s*v, and if there are
                                                   >> 705 //                   two intersections, s[0] < s[1]. May be negative.
                                                   >> 706 // Returns:
                                                   >> 707 //     The number of intersections. If 0, the trajectory misses. If 1, the 
                                                   >> 708 //     trajectory just grazes the surface.
                                                   >> 709 //
                                                   >> 710 // Solution:
                                                   >> 711 //     One needs to solve: ( (p.x + s*v.x)/dx )**2  + ( (p.y + s*v.y)/dy )**2 = 1
                                                   >> 712 //
                                                   >> 713 //     The solution is quadratic: a*s**2 + b*s + c = 0
                                                   >> 714 //
                                                   >> 715 //           a = (v.x/dx)**2 + (v.y/dy)**2
                                                   >> 716 //           b = 2*p.x*v.x/dx**2 + 2*p.y*v.y/dy**2
                                                   >> 717 //           c = (p.x/dx)**2 + (p.y/dy)**2 - 1
                                                   >> 718 //
                                                   >> 719 G4int G4EllipticalTube::IntersectXY( const G4ThreeVector &p,
                                                   >> 720                                      const G4ThreeVector &v,
                                                   >> 721                                            G4double s[2] ) const
                                                   >> 722 {
                                                   >> 723   G4double px = p.x(), py = p.y();
                                                   >> 724   G4double vx = v.x(), vy = v.y();
                                                   >> 725   
                                                   >> 726   G4double a = (vx/dx)*(vx/dx) + (vy/dy)*(vy/dy);
                                                   >> 727   G4double b = 2.0*( px*vx/dx/dx + py*vy/dy/dy );
                                                   >> 728   G4double c = (px/dx)*(px/dx) + (py/dy)*(py/dy) - 1.0;
                                                   >> 729   
                                                   >> 730   if (a < DBL_MIN) return 0;      // Trajectory parallel to z axis
                                                   >> 731   
                                                   >> 732   G4double radical = b*b - 4*a*c;
                                                   >> 733   
                                                   >> 734   if (radical < -DBL_MIN) return 0;    // No solution
                                                   >> 735   
                                                   >> 736   if (radical < DBL_MIN)
589   {                                               737   {
590     std::ostringstream message;                << 738     //
591     G4long oldprc = message.precision(16);     << 739     // Grazes surface
592     message << "Point p is outside (!?) of sol << 740     //
593             << "Position:\n"                   << 741     s[0] = -b/a/2.0;
594             << "   p.x() = "  << p.x()/mm << " << 742     return 1;
595             << "   p.y() = "  << p.y()/mm << " << 743   }
596             << "   p.z() = "  << p.z()/mm << " << 744   
597     message.precision(oldprc) ;                << 745   radical = std::sqrt(radical);
598     G4Exception("G4ElliptocalTube::DistanceToO << 746   
599                 JustWarning, message);         << 747   G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
600     DumpInfo();                                << 748   G4double sa = q/a;
601   }                                            << 749   G4double sb = c/q;    
602 #endif                                         << 750   if (sa < sb) { s[0] = sa; s[1] = sb; } else { s[0] = sb; s[1] = sa; }
603   // safety distance to Z-bases                << 751   return 2;
604   G4double distZ = fDz - std::abs(p.z());      << 
605                                                << 
606   // safety distance lateral surface           << 
607   G4double x = p.x() * fSx;                    << 
608   G4double y = p.y() * fSy;                    << 
609   G4double distR = fR - std::sqrt(x * x + y *  << 
610                                                << 
611   // return SafetyToOut                        << 
612   G4double dist = std::min(distZ, distR);      << 
613   return (dist < 0) ? 0 : dist;                << 
614 }                                                 752 }
615                                                   753 
616 ////////////////////////////////////////////// << 754 
617 //                                                755 //
618 // GetEntityType                                  756 // GetEntityType
619                                                << 757 //
620 G4GeometryType G4EllipticalTube::GetEntityType    758 G4GeometryType G4EllipticalTube::GetEntityType() const
621 {                                                 759 {
622   return {"G4EllipticalTube"};                 << 760   return G4String("G4EllipticalTube");
623 }                                                 761 }
624                                                   762 
625 ////////////////////////////////////////////// << 
626 //                                             << 
627 // Make a clone of the object                  << 
628                                                   763 
629 G4VSolid* G4EllipticalTube::Clone() const      << 
630 {                                              << 
631   return new G4EllipticalTube(*this);          << 
632 }                                              << 
633                                                << 
634 ////////////////////////////////////////////// << 
635 //                                                764 //
636 // Return volume                               << 765 // GetCubicVolume
637                                                << 766 //
638 G4double G4EllipticalTube::GetCubicVolume()       767 G4double G4EllipticalTube::GetCubicVolume()
639 {                                                 768 {
640   if (fCubicVolume == 0.)                      << 769   if(fCubicVolume != 0.) ;
641   {                                            << 770     else fCubicVolume = G4VSolid::GetCubicVolume(); 
642     fCubicVolume = twopi * fDx * fDy * fDz;    << 
643   }                                            << 
644   return fCubicVolume;                            771   return fCubicVolume;
645 }                                                 772 }
646                                                   773 
647 ////////////////////////////////////////////// << 
648 //                                             << 
649 // Return cached surface area                  << 
650                                                << 
651 G4double G4EllipticalTube::GetCachedSurfaceAre << 
652 {                                              << 
653   G4ThreadLocalStatic G4double cached_Dx = 0;  << 
654   G4ThreadLocalStatic G4double cached_Dy = 0;  << 
655   G4ThreadLocalStatic G4double cached_Dz = 0;  << 
656   G4ThreadLocalStatic G4double cached_area = 0 << 
657   if (cached_Dx != fDx || cached_Dy != fDy ||  << 
658   {                                            << 
659     cached_Dx = fDx;                           << 
660     cached_Dy = fDy;                           << 
661     cached_Dz = fDz;                           << 
662     cached_area = 2.*(pi*fDx*fDy + G4GeomTools << 
663   }                                            << 
664   return cached_area;                          << 
665 }                                              << 
666                                                   774 
667 ////////////////////////////////////////////// << 
668 //                                                775 //
669 // Return surface area                         << 776 // Stream object contents to an output stream
670                                                << 
671 G4double G4EllipticalTube::GetSurfaceArea()    << 
672 {                                              << 
673   if(fSurfaceArea == 0.)                       << 
674   {                                            << 
675     fSurfaceArea = GetCachedSurfaceArea();     << 
676   }                                            << 
677   return fSurfaceArea;                         << 
678 }                                              << 
679                                                << 
680 ////////////////////////////////////////////// << 
681 //                                                777 //
682 // Stream object contents to output stream     << 
683                                                << 
684 std::ostream& G4EllipticalTube::StreamInfo(std    778 std::ostream& G4EllipticalTube::StreamInfo(std::ostream& os) const
685 {                                                 779 {
686   G4long oldprc = os.precision(16);            << 
687   os << "-------------------------------------    780   os << "-----------------------------------------------------------\n"
688      << "    *** Dump for solid - " << GetName    781      << "    *** Dump for solid - " << GetName() << " ***\n"
689      << "    =================================    782      << "    ===================================================\n"
690      << " Solid type: G4EllipticalTube\n"         783      << " Solid type: G4EllipticalTube\n"
691      << " Parameters: \n"                         784      << " Parameters: \n"
692      << "    length Z: " << fDz/mm << " mm \n" << 785      << "    length Z: " << dz/mm << " mm \n"
693      << "    lateral surface equation: \n"     << 786      << "    surface equation in X and Y: \n"
694      << "       (X / " << fDx << ")^2 + (Y / " << 787      << "       (X / " << dx << ")^2 + (Y / " << dy << ")^2 = 1 \n"
695      << "-------------------------------------    788      << "-----------------------------------------------------------\n";
696   os.precision(oldprc);                        << 
697                                                   789 
698   return os;                                      790   return os;
699 }                                                 791 }
700                                                   792 
701                                                   793 
702 ////////////////////////////////////////////// << 
703 //                                                794 //
704 // Pick up a random point on the surface       << 795 // GetPointOnSurface
705                                                << 796 //
                                                   >> 797 // Randomly generates a point on the surface, 
                                                   >> 798 // with ~ uniform distribution across surface.
                                                   >> 799 //
706 G4ThreeVector G4EllipticalTube::GetPointOnSurf    800 G4ThreeVector G4EllipticalTube::GetPointOnSurface() const
707 {                                                 801 {
708   // Select surface (0 - base at -Z, 1 - base  << 802   G4double xRand, yRand, zRand, phi, cosphi, sinphi, zArea, cArea,p, chose;
709   //                                           << 
710   G4double sbase = pi * fDx * fDy;             << 
711   G4double ssurf = GetCachedSurfaceArea();     << 
712   G4double select = ssurf * G4UniformRand();   << 
713                                                   803 
714   G4int k = 0;                                 << 804   phi    = RandFlat::shoot(0., 2.*pi);
715   if (select > sbase) k = 1;                   << 805   cosphi = std::cos(phi);
716   if (select > 2. * sbase) k = 2;              << 806   sinphi = std::sin(phi);
717                                                << 807   
718   // Pick random point on selected surface (re << 808   // the ellipse perimeter from: "http://mathworld.wolfram.com/Ellipse.html"
719   //                                           << 809   //   m = (dx - dy)/(dx + dy);
720   G4ThreeVector p;                             << 810   //   k = 1.+1./4.*m*m+1./64.*sqr(m)*sqr(m)+1./256.*sqr(m)*sqr(m)*sqr(m);
721   switch (k) {                                 << 811   //   p = pi*(a+b)*k;
722     case 0: // base at -Z                      << 812 
723     {                                          << 813   // perimeter below from "http://www.efunda.com/math/areas/EllipseGen.cfm"
724       G4TwoVector rho = G4RandomPointInEllipse << 814 
725       p.set(rho.x(), rho.y(), -fDz);           << 815   p = 2.*pi*std::sqrt(0.5*(dx*dx+dy*dy));
726       break;                                   << 816 
727     }                                          << 817   cArea = 2.*dz*p;
728     case 1: // base at +Z                      << 818   zArea = pi*dx*dy;
729     {                                          << 819 
730       G4TwoVector rho = G4RandomPointInEllipse << 820   xRand = dx*cosphi;
731       p.set(rho.x(), rho.y(), fDz);            << 821   yRand = dy*sinphi;
732       break;                                   << 822   zRand = RandFlat::shoot(dz, -1.*dz);
733     }                                          << 823     
734     case 2: // lateral surface                 << 824   chose = RandFlat::shoot(0.,2.*zArea+cArea);
735     {                                          << 825   
736       G4TwoVector rho = G4RandomPointOnEllipse << 826   if( (chose>=0) && (chose < cArea) )
737       p.set(rho.x(), rho.y(), (2. * G4UniformR << 827   {
738       break;                                   << 828     return G4ThreeVector (xRand,yRand,zRand);
739     }                                          << 829   }
                                                   >> 830   else if( (chose >= cArea) && (chose < cArea + zArea) )
                                                   >> 831   {
                                                   >> 832     xRand = RandFlat::shoot(-1.*dx,dx);
                                                   >> 833     yRand = std::sqrt(1.-sqr(xRand/dx));
                                                   >> 834     yRand = RandFlat::shoot(-1.*yRand, yRand);
                                                   >> 835     return G4ThreeVector (xRand,yRand,dz); 
                                                   >> 836   }
                                                   >> 837   else
                                                   >> 838   { 
                                                   >> 839     xRand = RandFlat::shoot(-1.*dx,dx);
                                                   >> 840     yRand = std::sqrt(1.-sqr(xRand/dx));
                                                   >> 841     yRand = RandFlat::shoot(-1.*yRand, yRand);
                                                   >> 842     return G4ThreeVector (xRand,yRand,-1.*dz);
740   }                                               843   }
741   return p;                                    << 
742 }                                                 844 }
743                                                   845 
744                                                   846 
745 ////////////////////////////////////////////// << 
746 //                                                847 //
747 // CreatePolyhedron                               848 // CreatePolyhedron
748                                                << 849 //
749 G4Polyhedron* G4EllipticalTube::CreatePolyhedr    850 G4Polyhedron* G4EllipticalTube::CreatePolyhedron() const
750 {                                                 851 {
751   // create cylinder with radius=1...             852   // create cylinder with radius=1...
752   //                                              853   //
753   G4Polyhedron* eTube = new G4PolyhedronTube(0 << 854   G4Polyhedron* eTube = new G4PolyhedronTube(0.,1.,dz);
754                                                   855 
755   // apply non-uniform scaling...                 856   // apply non-uniform scaling...
756   //                                              857   //
757   eTube->Transform(G4Scale3D(fDx, fDy, 1.));   << 858   eTube->Transform(G4Scale3D(dx,dy,1.));
758   return eTube;                                << 859   return  eTube;
759 }                                                 860 }
760                                                   861 
761 ////////////////////////////////////////////// << 862 
762 //                                                863 //
763 // GetPolyhedron                                  864 // GetPolyhedron
764                                                << 865 //
765 G4Polyhedron* G4EllipticalTube::GetPolyhedron     866 G4Polyhedron* G4EllipticalTube::GetPolyhedron () const
766 {                                                 867 {
767   if (fpPolyhedron == nullptr ||               << 868   if (!fpPolyhedron ||
768       fRebuildPolyhedron ||                    << 
769       fpPolyhedron->GetNumberOfRotationStepsAt    869       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
770       fpPolyhedron->GetNumberOfRotationSteps()    870       fpPolyhedron->GetNumberOfRotationSteps())
771   {                                            << 871     {
772     G4AutoLock l(&polyhedronMutex);            << 872       delete fpPolyhedron;
773     delete fpPolyhedron;                       << 873       fpPolyhedron = CreatePolyhedron();
774     fpPolyhedron = CreatePolyhedron();         << 874     }
775     fRebuildPolyhedron = false;                << 
776     l.unlock();                                << 
777   }                                            << 
778   return fpPolyhedron;                            875   return fpPolyhedron;
779 }                                                 876 }
780                                                   877 
781 ////////////////////////////////////////////// << 878 
782 //                                                879 //
783 // DescribeYourselfTo                             880 // DescribeYourselfTo
784                                                << 881 //
785 void G4EllipticalTube::DescribeYourselfTo( G4V    882 void G4EllipticalTube::DescribeYourselfTo( G4VGraphicsScene& scene ) const
786 {                                                 883 {
787   scene.AddSolid (*this);                         884   scene.AddSolid (*this);
788 }                                                 885 }
789                                                   886 
790 ////////////////////////////////////////////// << 887 
791 //                                                888 //
792 // GetExtent                                      889 // GetExtent
793                                                << 890 //
794 G4VisExtent G4EllipticalTube::GetExtent() cons    891 G4VisExtent G4EllipticalTube::GetExtent() const
795 {                                                 892 {
796   return { -fDx, fDx, -fDy, fDy, -fDz, fDz };  << 893   return G4VisExtent( -dx, dx, -dy, dy, -dz, dz );
797 }                                                 894 }
798                                                << 
799 #endif // !defined(G4GEOM_USE_UELLIPTICALTUBE) << 
800                                                   895