Geant4 Cross Reference

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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4Paraboloid.cc 83572 2014-09-01 15:23:27Z gcosmo $
                                                   >>  27 //
                                                   >>  28 // class G4Paraboloid
                                                   >>  29 //
 26 // Implementation for G4Paraboloid class           30 // Implementation for G4Paraboloid class
 27 //                                                 31 //
 28 // Author : Lukas Lindroos (CERN), July 2007       32 // Author : Lukas Lindroos (CERN), July 2007
 29 // Revised: Tatiana Nikitina (CERN)                33 // Revised: Tatiana Nikitina (CERN)
 30 // -------------------------------------------     34 // --------------------------------------------------------------------
 31                                                    35 
 32 #include "globals.hh"                              36 #include "globals.hh"
 33                                                    37 
 34 #include "G4Paraboloid.hh"                         38 #include "G4Paraboloid.hh"
 35                                                    39 
 36 #if !(defined(G4GEOM_USE_UPARABOLOID) && defin << 
 37                                                << 
 38 #include "G4VoxelLimits.hh"                        40 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"                    41 #include "G4AffineTransform.hh"
 40 #include "G4BoundingEnvelope.hh"               << 
 41                                                    42 
 42 #include "meshdefs.hh"                             43 #include "meshdefs.hh"
 43                                                    44 
 44 #include "Randomize.hh"                            45 #include "Randomize.hh"
 45                                                    46 
 46 #include "G4VGraphicsScene.hh"                     47 #include "G4VGraphicsScene.hh"
 47 #include "G4VisExtent.hh"                          48 #include "G4VisExtent.hh"
 48                                                    49 
 49 #include "G4AutoLock.hh"                           50 #include "G4AutoLock.hh"
 50                                                    51 
 51 namespace                                          52 namespace
 52 {                                                  53 {
 53   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE     54   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
 54 }                                                  55 }
 55                                                    56 
 56 using namespace CLHEP;                             57 using namespace CLHEP;
 57                                                    58 
 58 //////////////////////////////////////////////     59 ///////////////////////////////////////////////////////////////////////////////
 59 //                                                 60 //
 60 // constructor - check parameters                  61 // constructor - check parameters
 61 //                                             <<  62 
 62 G4Paraboloid::G4Paraboloid(const G4String& pNa     63 G4Paraboloid::G4Paraboloid(const G4String& pName,
 63                                  G4double pDz,     64                                  G4double pDz,
 64                                  G4double pR1,     65                                  G4double pR1,
 65                                  G4double pR2)     66                                  G4double pR2)
 66  : G4VSolid(pName)                             <<  67  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
                                                   >>  68    fSurfaceArea(0.), fCubicVolume(0.) 
                                                   >>  69 
 67 {                                                  70 {
 68   if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.     71   if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.) )
 69   {                                                72   {
 70     std::ostringstream message;                    73     std::ostringstream message;
 71     message << "Invalid dimensions. Negative I     74     message << "Invalid dimensions. Negative Input Values or R1>=R2 - "
 72             << GetName();                          75             << GetName();
 73     G4Exception("G4Paraboloid::G4Paraboloid()"     76     G4Exception("G4Paraboloid::G4Paraboloid()", "GeomSolids0002", 
 74                 FatalErrorInArgument, message,     77                 FatalErrorInArgument, message,
 75                 "Z half-length must be larger      78                 "Z half-length must be larger than zero or R1>=R2.");
 76   }                                                79   }
 77                                                    80 
 78   r1 = pR1;                                        81   r1 = pR1;
 79   r2 = pR2;                                        82   r2 = pR2;
 80   dz = pDz;                                        83   dz = pDz;
 81                                                    84 
 82   // r1^2 = k1 * (-dz) + k2                        85   // r1^2 = k1 * (-dz) + k2
 83   // r2^2 = k1 * ( dz) + k2                        86   // r2^2 = k1 * ( dz) + k2
 84   // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 +      87   // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2
 85   // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) =>     88   // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz
 86                                                    89 
 87   k1 = (r2 * r2 - r1 * r1) / 2 / dz;               90   k1 = (r2 * r2 - r1 * r1) / 2 / dz;
 88   k2 = (r2 * r2 + r1 * r1) / 2;                    91   k2 = (r2 * r2 + r1 * r1) / 2;
 89 }                                                  92 }
 90                                                    93 
 91 //////////////////////////////////////////////     94 ///////////////////////////////////////////////////////////////////////////////
 92 //                                                 95 //
 93 // Fake default constructor - sets only member     96 // Fake default constructor - sets only member data and allocates memory
 94 //                            for usage restri     97 //                            for usage restricted to object persistency.
 95 //                                                 98 //
 96 G4Paraboloid::G4Paraboloid( __void__& a )          99 G4Paraboloid::G4Paraboloid( __void__& a )
 97   : G4VSolid(a), dz(0.), r1(0.), r2(0.), k1(0. << 100   : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
                                                   >> 101     fSurfaceArea(0.), fCubicVolume(0.),
                                                   >> 102     dz(0.), r1(0.), r2(0.), k1(0.), k2(0.)
 98 {                                                 103 {
 99 }                                                 104 }
100                                                   105 
101 //////////////////////////////////////////////    106 ///////////////////////////////////////////////////////////////////////////////
102 //                                                107 //
103 // Destructor                                     108 // Destructor
104 //                                             << 109 
105 G4Paraboloid::~G4Paraboloid()                     110 G4Paraboloid::~G4Paraboloid()
106 {                                                 111 {
107   delete fpPolyhedron; fpPolyhedron = nullptr; << 112   delete fpPolyhedron; fpPolyhedron = 0;
108 }                                                 113 }
109                                                   114 
110 //////////////////////////////////////////////    115 ///////////////////////////////////////////////////////////////////////////////
111 //                                                116 //
112 // Copy constructor                               117 // Copy constructor
113 //                                             << 118 
114 G4Paraboloid::G4Paraboloid(const G4Paraboloid&    119 G4Paraboloid::G4Paraboloid(const G4Paraboloid& rhs)
115   : G4VSolid(rhs),                             << 120   : G4VSolid(rhs), fRebuildPolyhedron(false), fpPolyhedron(0),
116     fSurfaceArea(rhs.fSurfaceArea), fCubicVolu    121     fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume),
117     dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs    122     dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs.k1), k2(rhs.k2)
118 {                                                 123 {
119 }                                                 124 }
120                                                   125 
                                                   >> 126 
121 //////////////////////////////////////////////    127 ///////////////////////////////////////////////////////////////////////////////
122 //                                                128 //
123 // Assignment operator                            129 // Assignment operator
124 //                                             << 130 
125 G4Paraboloid& G4Paraboloid::operator = (const     131 G4Paraboloid& G4Paraboloid::operator = (const G4Paraboloid& rhs) 
126 {                                                 132 {
127    // Check assignment to self                    133    // Check assignment to self
128    //                                             134    //
129    if (this == &rhs)  { return *this; }           135    if (this == &rhs)  { return *this; }
130                                                   136 
131    // Copy base class data                        137    // Copy base class data
132    //                                             138    //
133    G4VSolid::operator=(rhs);                      139    G4VSolid::operator=(rhs);
134                                                   140 
135    // Copy data                                   141    // Copy data
136    //                                             142    //
137    fSurfaceArea = rhs.fSurfaceArea; fCubicVolu    143    fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
138    dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 =    144    dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = rhs.k1; k2 = rhs.k2;
139    fRebuildPolyhedron = false;                    145    fRebuildPolyhedron = false;
140    delete fpPolyhedron; fpPolyhedron = nullptr << 146    delete fpPolyhedron; fpPolyhedron = 0;
141                                                   147 
142    return *this;                                  148    return *this;
143 }                                                 149 }
144                                                   150 
145 ////////////////////////////////////////////// << 151 /////////////////////////////////////////////////////////////////////////
146 //                                             << 
147 // Get bounding box                            << 
148 //                                                152 //
149 void G4Paraboloid::BoundingLimits(G4ThreeVecto << 153 // Dispatch to parameterisation for replication mechanism dimension
150                                   G4ThreeVecto << 154 // computation & modification.
151 {                                              << 155 
152   pMin.set(-r2,-r2,-dz);                       << 156 //void ComputeDimensions(       G4VPVParamerisation p,
153   pMax.set( r2, r2, dz);                       << 157 //                        const G4Int               n,
                                                   >> 158 //                        const G4VPhysicalVolume*  pRep )
                                                   >> 159 //{
                                                   >> 160 //  p->ComputeDimensions(*this,n,pRep) ;
                                                   >> 161 //}
154                                                   162 
155   // Check correctness of the bounding box     << 
156   //                                           << 
157   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
158   {                                            << 
159     std::ostringstream message;                << 
160     message << "Bad bounding box (min >= max)  << 
161             << GetName() << " !"               << 
162             << "\npMin = " << pMin             << 
163             << "\npMax = " << pMax;            << 
164     G4Exception("G4Paraboloid::BoundingLimits( << 
165                 JustWarning, message);         << 
166     DumpInfo();                                << 
167   }                                            << 
168 }                                              << 
169                                                   163 
170 //////////////////////////////////////////////    164 ///////////////////////////////////////////////////////////////////////////////
171 //                                                165 //
172 // Calculate extent under transform and specif    166 // Calculate extent under transform and specified limit
173 //                                             << 167 
174 G4bool                                            168 G4bool
175 G4Paraboloid::CalculateExtent(const EAxis pAxi    169 G4Paraboloid::CalculateExtent(const EAxis pAxis,
176                               const G4VoxelLim << 170                              const G4VoxelLimits& pVoxelLimit,
177                               const G4AffineTr << 171                              const G4AffineTransform& pTransform,
178                                     G4double&  << 172                                    G4double& pMin, G4double& pMax) const
179 {                                              << 173 {
180   G4ThreeVector bmin, bmax;                    << 174   G4double xMin = -r2 + pTransform.NetTranslation().x(),
181                                                << 175            xMax = r2 + pTransform.NetTranslation().x(),
182   // Get bounding box                          << 176            yMin = -r2 + pTransform.NetTranslation().y(),
183   BoundingLimits(bmin,bmax);                   << 177            yMax = r2 + pTransform.NetTranslation().y(),
184                                                << 178            zMin = -dz + pTransform.NetTranslation().z(),
185   // Find extent                               << 179            zMax = dz + pTransform.NetTranslation().z();
186   G4BoundingEnvelope bbox(bmin,bmax);          << 180 
187   return bbox.CalculateExtent(pAxis,pVoxelLimi << 181   if(!pTransform.IsRotated()
                                                   >> 182   || pTransform.NetRotation()(G4ThreeVector(0, 0, 1)) == G4ThreeVector(0, 0, 1))
                                                   >> 183   {
                                                   >> 184     if(pVoxelLimit.IsXLimited())
                                                   >> 185     {
                                                   >> 186       if(pVoxelLimit.GetMaxXExtent() < xMin - 0.5 * kCarTolerance
                                                   >> 187       || pVoxelLimit.GetMinXExtent() > xMax + 0.5 * kCarTolerance)
                                                   >> 188       {
                                                   >> 189         return false;
                                                   >> 190       }
                                                   >> 191       else
                                                   >> 192       {
                                                   >> 193         if(pVoxelLimit.GetMinXExtent() > xMin)
                                                   >> 194         {
                                                   >> 195           xMin = pVoxelLimit.GetMinXExtent();
                                                   >> 196         }
                                                   >> 197         if(pVoxelLimit.GetMaxXExtent() < xMax)
                                                   >> 198         {
                                                   >> 199           xMax = pVoxelLimit.GetMaxXExtent();
                                                   >> 200         }
                                                   >> 201       }
                                                   >> 202     }
                                                   >> 203     if(pVoxelLimit.IsYLimited())
                                                   >> 204     {
                                                   >> 205       if(pVoxelLimit.GetMaxYExtent() < yMin - 0.5 * kCarTolerance
                                                   >> 206       || pVoxelLimit.GetMinYExtent() > yMax + 0.5 * kCarTolerance)
                                                   >> 207       {
                                                   >> 208         return false;
                                                   >> 209       }
                                                   >> 210       else
                                                   >> 211       {
                                                   >> 212         if(pVoxelLimit.GetMinYExtent() > yMin)
                                                   >> 213         {
                                                   >> 214           yMin = pVoxelLimit.GetMinYExtent();
                                                   >> 215         }
                                                   >> 216         if(pVoxelLimit.GetMaxYExtent() < yMax)
                                                   >> 217         {
                                                   >> 218           yMax = pVoxelLimit.GetMaxYExtent();
                                                   >> 219         }
                                                   >> 220       }
                                                   >> 221     }
                                                   >> 222     if(pVoxelLimit.IsZLimited())
                                                   >> 223     {
                                                   >> 224       if(pVoxelLimit.GetMaxZExtent() < zMin - 0.5 * kCarTolerance
                                                   >> 225       || pVoxelLimit.GetMinZExtent() > zMax + 0.5 * kCarTolerance)
                                                   >> 226       {
                                                   >> 227         return false;
                                                   >> 228       }
                                                   >> 229       else
                                                   >> 230       {
                                                   >> 231         if(pVoxelLimit.GetMinZExtent() > zMin)
                                                   >> 232         {
                                                   >> 233           zMin = pVoxelLimit.GetMinZExtent();
                                                   >> 234         }
                                                   >> 235         if(pVoxelLimit.GetMaxZExtent() < zMax)
                                                   >> 236         {
                                                   >> 237           zMax = pVoxelLimit.GetMaxZExtent();
                                                   >> 238         }
                                                   >> 239       }
                                                   >> 240     }
                                                   >> 241     switch(pAxis)
                                                   >> 242     {
                                                   >> 243       case kXAxis:
                                                   >> 244         pMin = xMin;
                                                   >> 245         pMax = xMax;
                                                   >> 246         break;
                                                   >> 247       case kYAxis:
                                                   >> 248         pMin = yMin;
                                                   >> 249         pMax = yMax;
                                                   >> 250         break;
                                                   >> 251       case kZAxis:
                                                   >> 252         pMin = zMin;
                                                   >> 253         pMax = zMax;
                                                   >> 254         break;
                                                   >> 255       default:
                                                   >> 256         pMin = 0;
                                                   >> 257         pMax = 0;
                                                   >> 258         return false;
                                                   >> 259     }
                                                   >> 260   }
                                                   >> 261   else
                                                   >> 262   {
                                                   >> 263     G4bool existsAfterClip=true;
                                                   >> 264 
                                                   >> 265     // Calculate rotated vertex coordinates
                                                   >> 266 
                                                   >> 267     G4int noPolygonVertices=0;
                                                   >> 268     G4ThreeVectorList* vertices
                                                   >> 269       = CreateRotatedVertices(pTransform,noPolygonVertices);
                                                   >> 270 
                                                   >> 271     if(pAxis == kXAxis || pAxis == kYAxis || pAxis == kZAxis)
                                                   >> 272     {
                                                   >> 273 
                                                   >> 274       pMin =  kInfinity;
                                                   >> 275       pMax = -kInfinity;
                                                   >> 276 
                                                   >> 277       for(G4ThreeVectorList::iterator it = vertices->begin();
                                                   >> 278           it < vertices->end(); it++)
                                                   >> 279       {
                                                   >> 280         if(pMin > (*it)[pAxis]) pMin = (*it)[pAxis];
                                                   >> 281         if((*it)[pAxis] < pVoxelLimit.GetMinExtent(pAxis))
                                                   >> 282         {
                                                   >> 283           pMin = pVoxelLimit.GetMinExtent(pAxis);
                                                   >> 284         }
                                                   >> 285         if(pMax < (*it)[pAxis])
                                                   >> 286         {
                                                   >> 287           pMax = (*it)[pAxis];
                                                   >> 288         }
                                                   >> 289         if((*it)[pAxis] > pVoxelLimit.GetMaxExtent(pAxis))
                                                   >> 290         {
                                                   >> 291           pMax = pVoxelLimit.GetMaxExtent(pAxis);
                                                   >> 292         }
                                                   >> 293       }
                                                   >> 294 
                                                   >> 295       if(pMin > pVoxelLimit.GetMaxExtent(pAxis)
                                                   >> 296       || pMax < pVoxelLimit.GetMinExtent(pAxis)) { existsAfterClip = false; }
                                                   >> 297     }
                                                   >> 298     else
                                                   >> 299     {
                                                   >> 300       pMin = 0;
                                                   >> 301       pMax = 0;
                                                   >> 302       existsAfterClip = false;
                                                   >> 303     }
                                                   >> 304     delete vertices;
                                                   >> 305     return existsAfterClip;
                                                   >> 306   }
                                                   >> 307   return true;
188 }                                                 308 }
189                                                   309 
190 //////////////////////////////////////////////    310 ///////////////////////////////////////////////////////////////////////////////
191 //                                                311 //
192 // Return whether point inside/outside/on surf    312 // Return whether point inside/outside/on surface
193 //                                             << 313 
194 EInside G4Paraboloid::Inside(const G4ThreeVect    314 EInside G4Paraboloid::Inside(const G4ThreeVector& p) const
195 {                                                 315 {
196   // First check is  the point is above or bel    316   // First check is  the point is above or below the solid.
197   //                                              317   //
198   if(std::fabs(p.z()) > dz + 0.5 * kCarToleran    318   if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; }
199                                                   319 
200   G4double rho2 = p.perp2(),                      320   G4double rho2 = p.perp2(),
201            rhoSurfTimesTol2  = (k1 * p.z() + k    321            rhoSurfTimesTol2  = (k1 * p.z() + k2) * sqr(kCarTolerance),
202            A = rho2 - ((k1 *p.z() + k2) + 0.25    322            A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance);
203                                                   323  
204   if(A < 0 && sqr(A) > rhoSurfTimesTol2)          324   if(A < 0 && sqr(A) > rhoSurfTimesTol2)
205   {                                               325   {
206     // Actually checking rho < radius of parab    326     // Actually checking rho < radius of paraboloid at z = p.z().
207     // We're either inside or in lower/upper c    327     // We're either inside or in lower/upper cutoff area.
208                                                   328    
209     if(std::fabs(p.z()) > dz - 0.5 * kCarToler    329     if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
210     {                                             330     {
211       // We're in the upper/lower cutoff area,    331       // We're in the upper/lower cutoff area, sides have a paraboloid shape
212       // maybe further checks should be made t    332       // maybe further checks should be made to make these nicer
213                                                   333 
214       return kSurface;                            334       return kSurface;
215     }                                             335     }
216     else                                          336     else
217     {                                             337     {
218       return kInside;                             338       return kInside;
219     }                                             339     }
220   }                                               340   }
221   else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)    341   else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
222   {                                               342   {
223     // We're in the parabolic surface.            343     // We're in the parabolic surface.
224                                                   344 
225     return kSurface;                              345     return kSurface;
226   }                                               346   }
227   else                                            347   else
228   {                                               348   {
229     return kOutside;                              349     return kOutside;
230   }                                               350   }
231 }                                                 351 }
232                                                   352 
233 //////////////////////////////////////////////    353 ///////////////////////////////////////////////////////////////////////////////
234 //                                                354 //
235 // SurfaceNormal                               << 355 
236 //                                             << 
237 G4ThreeVector G4Paraboloid::SurfaceNormal( con    356 G4ThreeVector G4Paraboloid::SurfaceNormal( const G4ThreeVector& p) const
238 {                                                 357 {
239   G4ThreeVector n(0, 0, 0);                       358   G4ThreeVector n(0, 0, 0);
240   if(std::fabs(p.z()) > dz + 0.5 * kCarToleran    359   if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance)
241   {                                               360   {
242     // If above or below just return normal ve    361     // If above or below just return normal vector for the cutoff plane.
243                                                   362 
244     n = G4ThreeVector(0, 0, p.z()/std::fabs(p.    363     n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z()));
245   }                                               364   }
246   else if(std::fabs(p.z()) > dz - 0.5 * kCarTo    365   else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
247   {                                               366   {
248     // This means we're somewhere in the plane    367     // This means we're somewhere in the plane z = dz or z = -dz.
249     // (As far as the program is concerned any    368     // (As far as the program is concerned anyway.
250                                                   369 
251     if(p.z() < 0) // Are we in upper or lower     370     if(p.z() < 0) // Are we in upper or lower plane?
252     {                                             371     {
253       if(p.perp2() > sqr(r1 + 0.5 * kCarTolera    372       if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance))
254       {                                           373       {
255         n = G4ThreeVector(p.x(), p.y(), -k1 /     374         n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit();
256       }                                           375       }
257       else if(r1 < 0.5 * kCarTolerance            376       else if(r1 < 0.5 * kCarTolerance
258            || p.perp2() > sqr(r1 - 0.5 * kCarT    377            || p.perp2() > sqr(r1 - 0.5 * kCarTolerance))
259       {                                           378       {
260         n = G4ThreeVector(p.x(), p.y(), 0.).un    379         n = G4ThreeVector(p.x(), p.y(), 0.).unit()
261           + G4ThreeVector(0., 0., -1.).unit();    380           + G4ThreeVector(0., 0., -1.).unit();
262         n = n.unit();                             381         n = n.unit();
263       }                                           382       }
264       else                                        383       else
265       {                                           384       {
266         n = G4ThreeVector(0., 0., -1.);           385         n = G4ThreeVector(0., 0., -1.);
267       }                                           386       }
268     }                                             387     }
269     else                                          388     else
270     {                                             389     {
271       if(p.perp2() > sqr(r2 + 0.5 * kCarTolera    390       if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance))
272       {                                           391       {
273         n = G4ThreeVector(p.x(), p.y(), 0.).un    392         n = G4ThreeVector(p.x(), p.y(), 0.).unit();
274       }                                           393       }
275       else if(r2 < 0.5 * kCarTolerance            394       else if(r2 < 0.5 * kCarTolerance
276            || p.perp2() > sqr(r2 - 0.5 * kCarT    395            || p.perp2() > sqr(r2 - 0.5 * kCarTolerance))
277       {                                           396       {
278         n = G4ThreeVector(p.x(), p.y(), 0.).un    397         n = G4ThreeVector(p.x(), p.y(), 0.).unit()
279           + G4ThreeVector(0., 0., 1.).unit();     398           + G4ThreeVector(0., 0., 1.).unit();
280         n = n.unit();                             399         n = n.unit();
281       }                                           400       }
282       else                                        401       else
283       {                                           402       {
284         n = G4ThreeVector(0., 0., 1.);            403         n = G4ThreeVector(0., 0., 1.);
285       }                                           404       }
286     }                                             405     }
287   }                                               406   }
288   else                                            407   else
289   {                                               408   {
290     G4double rho2 = p.perp2();                    409     G4double rho2 = p.perp2();
291     G4double rhoSurfTimesTol2  = (k1 * p.z() +    410     G4double rhoSurfTimesTol2  = (k1 * p.z() + k2) * sqr(kCarTolerance);
292     G4double A = rho2 - ((k1 *p.z() + k2)         411     G4double A = rho2 - ((k1 *p.z() + k2)
293                + 0.25 * kCarTolerance * kCarTo    412                + 0.25 * kCarTolerance * kCarTolerance);
294                                                   413 
295     if(A < 0 && sqr(A) > rhoSurfTimesTol2)        414     if(A < 0 && sqr(A) > rhoSurfTimesTol2)
296     {                                             415     {
297       // Actually checking rho < radius of par    416       // Actually checking rho < radius of paraboloid at z = p.z().
298       // We're inside.                            417       // We're inside.
299                                                   418 
300       if(p.mag2() != 0) { n = p.unit(); }         419       if(p.mag2() != 0) { n = p.unit(); }
301     }                                             420     }
302     else if(A <= 0 || sqr(A) < rhoSurfTimesTol    421     else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
303     {                                             422     {
304       // We're in the parabolic surface.          423       // We're in the parabolic surface.
305                                                   424 
306       n = G4ThreeVector(p.x(), p.y(), - k1 / 2    425       n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
307     }                                             426     }
308     else                                          427     else
309     {                                             428     {
310       n = G4ThreeVector(p.x(), p.y(), - k1 / 2    429       n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
311     }                                             430     }
312   }                                               431   }
313                                                   432 
314   if(n.mag2() == 0)                               433   if(n.mag2() == 0)
315   {                                               434   {
316     std::ostringstream message;                   435     std::ostringstream message;
317     message << "No normal defined for this poi    436     message << "No normal defined for this point p." << G4endl
318             << "          p = " << 1 / mm * p     437             << "          p = " << 1 / mm * p << " mm";
319     G4Exception("G4Paraboloid::SurfaceNormal(p    438     G4Exception("G4Paraboloid::SurfaceNormal(p)", "GeomSolids1002",
320                 JustWarning, message);            439                 JustWarning, message);
321   }                                               440   }
322   return n;                                       441   return n;
323 }                                                 442 }
324                                                   443 
325 //////////////////////////////////////////////    444 ///////////////////////////////////////////////////////////////////////////////
326 //                                                445 //
327 // Calculate distance to shape from outside, a    446 // Calculate distance to shape from outside, along normalised vector
328 // - return kInfinity if no intersection          447 // - return kInfinity if no intersection
329 //                                                448 //
330 //                                             << 449 
331 G4double G4Paraboloid::DistanceToIn( const G4T    450 G4double G4Paraboloid::DistanceToIn( const G4ThreeVector& p,
332                                      const G4T << 451                                     const G4ThreeVector& v  ) const
333 {                                                 452 {
334   G4double rho2 = p.perp2(), paraRho2 = std::f    453   G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
335   G4double tol2 = kCarTolerance*kCarTolerance;    454   G4double tol2 = kCarTolerance*kCarTolerance;
336   G4double tolh = 0.5*kCarTolerance;              455   G4double tolh = 0.5*kCarTolerance;
337                                                   456 
338   if((r2 != 0.0) && p.z() > - tolh + dz)       << 457   if(r2 && p.z() > - tolh + dz) 
339   {                                               458   {
340     // If the points is above check for inters    459     // If the points is above check for intersection with upper edge.
341                                                   460 
342     if(v.z() < 0)                                 461     if(v.z() < 0)
343     {                                             462     {
344       G4double intersection = (dz - p.z()) / v    463       G4double intersection = (dz - p.z()) / v.z(); // With plane z = dz.
345       if(sqr(p.x() + v.x()*intersection)          464       if(sqr(p.x() + v.x()*intersection)
346        + sqr(p.y() + v.y()*intersection) < sqr    465        + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance))
347       {                                           466       {
348         if(p.z() < tolh + dz)                     467         if(p.z() < tolh + dz)
349           { return 0; }                           468           { return 0; }
350         else                                      469         else
351           { return intersection; }                470           { return intersection; }
352       }                                           471       }
353     }                                             472     }
354     else  // Direction away, no possibility of    473     else  // Direction away, no possibility of intersection
355     {                                             474     {
356       return kInfinity;                           475       return kInfinity;
357     }                                             476     }
358   }                                               477   }
359   else if((r1 != 0.0) && p.z() < tolh - dz)    << 478   else if(r1 && p.z() < tolh - dz)
360   {                                               479   {
361     // If the points is belove check for inter    480     // If the points is belove check for intersection with lower edge.
362                                                   481 
363     if(v.z() > 0)                                 482     if(v.z() > 0)
364     {                                             483     {
365       G4double intersection = (-dz - p.z()) /     484       G4double intersection = (-dz - p.z()) / v.z(); // With plane z = -dz.
366       if(sqr(p.x() + v.x()*intersection)          485       if(sqr(p.x() + v.x()*intersection)
367        + sqr(p.y() + v.y()*intersection) < sqr    486        + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance))
368       {                                           487       {
369         if(p.z() > -tolh - dz)                    488         if(p.z() > -tolh - dz)
370         {                                         489         {
371           return 0;                               490           return 0;
372         }                                         491         }
373         else                                      492         else
374         {                                         493         {
375           return intersection;                    494           return intersection;
376         }                                         495         }
377       }                                           496       }
378     }                                             497     }
379     else  // Direction away, no possibility of    498     else  // Direction away, no possibility of intersection
380     {                                             499     {
381       return kInfinity;                           500       return kInfinity;
382     }                                             501     }
383   }                                               502   }
384                                                   503 
385   G4double A = k1 / 2 * v.z() - p.x() * v.x()     504   G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(),
386            vRho2 = v.perp2(), intersection,       505            vRho2 = v.perp2(), intersection,
387            B = (k1 * p.z() + k2 - rho2) * vRho    506            B = (k1 * p.z() + k2 - rho2) * vRho2;
388                                                   507 
389   if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRh    508   if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) )
390     || (p.z() < - dz+kCarTolerance)               509     || (p.z() < - dz+kCarTolerance)
391     || (p.z() > dz-kCarTolerance) ) // Make su    510     || (p.z() > dz-kCarTolerance) ) // Make sure it's safely outside.
392   {                                               511   {
393     // Is there a problem with squaring rho tw    512     // Is there a problem with squaring rho twice?
394                                                   513 
395     if(vRho2<tol2) // Needs to be treated sepe    514     if(vRho2<tol2) // Needs to be treated seperately.
396     {                                             515     {
397       intersection = ((rho2 - k2)/k1 - p.z())/    516       intersection = ((rho2 - k2)/k1 - p.z())/v.z();
398       if(intersection < 0) { return kInfinity;    517       if(intersection < 0) { return kInfinity; }
399       else if(std::fabs(p.z() + v.z() * inters    518       else if(std::fabs(p.z() + v.z() * intersection) <= dz)
400       {                                           519       {
401         return intersection;                      520         return intersection;
402       }                                           521       }
403       else                                        522       else
404       {                                           523       {
405         return kInfinity;                         524         return kInfinity;
406       }                                           525       }
407     }                                             526     }
408     else if(A*A + B < 0) // No real intersecti    527     else if(A*A + B < 0) // No real intersections.
409     {                                             528     {
410       return kInfinity;                           529       return kInfinity;
411     }                                             530     }
412     else                                          531     else
413     {                                             532     {
414       intersection = (A - std::sqrt(B + sqr(A)    533       intersection = (A - std::sqrt(B + sqr(A))) / vRho2;
415       if(intersection < 0)                        534       if(intersection < 0)
416       {                                           535       {
417         return kInfinity;                         536         return kInfinity;
418       }                                           537       }
419       else if(std::fabs(p.z() + intersection *    538       else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh)
420       {                                           539       {
421         return intersection;                      540         return intersection;
422       }                                           541       }
423       else                                        542       else
424       {                                           543       {
425         return kInfinity;                         544         return kInfinity;
426       }                                           545       }
427     }                                             546     }
428   }                                               547   }
429   else if(sqr(rho2 - paraRho2 - .25 * tol2) <=    548   else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2)
430   {                                               549   {
431     // If this is true we're somewhere in the     550     // If this is true we're somewhere in the border.
432                                                   551 
433     G4ThreeVector normal(p.x(), p.y(), -k1/2);    552     G4ThreeVector normal(p.x(), p.y(), -k1/2);
434     if(normal.dot(v) <= 0)                        553     if(normal.dot(v) <= 0)
435       { return 0; }                               554       { return 0; }
436     else                                          555     else
437       { return kInfinity; }                       556       { return kInfinity; }
438   }                                               557   }
439   else                                            558   else
440   {                                               559   {
441     std::ostringstream message;                   560     std::ostringstream message;
442     if(Inside(p) == kInside)                      561     if(Inside(p) == kInside)
443     {                                             562     {
444       message << "Point p is inside! - " << Ge    563       message << "Point p is inside! - " << GetName() << G4endl;
445     }                                             564     }
446     else                                          565     else
447     {                                             566     {
448       message << "Likely a problem in this fun    567       message << "Likely a problem in this function, for solid: " << GetName()
449               << G4endl;                          568               << G4endl;
450     }                                             569     }
451     message << "          p = " << p * (1/mm)     570     message << "          p = " << p * (1/mm) << " mm" << G4endl
452             << "          v = " << v * (1/mm)     571             << "          v = " << v * (1/mm) << " mm";
453     G4Exception("G4Paraboloid::DistanceToIn(p,    572     G4Exception("G4Paraboloid::DistanceToIn(p,v)", "GeomSolids1002",
454                 JustWarning, message);            573                 JustWarning, message);
455     return 0;                                     574     return 0;
456   }                                               575   }
457 }                                                 576 }
458                                                   577 
459 //////////////////////////////////////////////    578 ///////////////////////////////////////////////////////////////////////////////
460 //                                                579 //
461 // Calculate distance (<= actual) to closest s    580 // Calculate distance (<= actual) to closest surface of shape from outside
462 // - Return zero if point inside               << 581 // - Return 0 if point inside
463 //                                             << 582 
464 G4double G4Paraboloid::DistanceToIn(const G4Th    583 G4double G4Paraboloid::DistanceToIn(const G4ThreeVector& p) const
465 {                                                 584 {
466   G4double safz = -dz+std::fabs(p.z());           585   G4double safz = -dz+std::fabs(p.z());
467   if(safz<0.) { safz=0.; }                     << 586   if(safz<0) { safz=0; }
468   G4double safr = kInfinity;                      587   G4double safr = kInfinity;
469                                                   588 
470   G4double rho = p.x()*p.x()+p.y()*p.y();         589   G4double rho = p.x()*p.x()+p.y()*p.y();
471   G4double paraRho = (p.z()-k2)/k1;               590   G4double paraRho = (p.z()-k2)/k1;
472   G4double sqrho = std::sqrt(rho);                591   G4double sqrho = std::sqrt(rho);
473                                                   592 
474   if(paraRho<0.)                               << 593   if(paraRho<0)
475   {                                               594   {
476     safr=sqrho-r2;                                595     safr=sqrho-r2;
477     if(safr>safz) { safz=safr; }                  596     if(safr>safz) { safz=safr; }
478     return safz;                                  597     return safz;
479   }                                               598   }
480                                                   599 
481   G4double sqprho = std::sqrt(paraRho);           600   G4double sqprho = std::sqrt(paraRho);
482   G4double dRho = sqrho-sqprho;                   601   G4double dRho = sqrho-sqprho;
483   if(dRho<0.) { return safz; }                 << 602   if(dRho<0) { return safz; }
484                                                   603 
485   G4double talf = -2.*k1*sqprho;                  604   G4double talf = -2.*k1*sqprho;
486   G4double tmp  = 1+talf*talf;                    605   G4double tmp  = 1+talf*talf;
487   if(tmp<0.) { return safz; }                     606   if(tmp<0.) { return safz; }
488                                                   607 
489   G4double salf = talf/std::sqrt(tmp);            608   G4double salf = talf/std::sqrt(tmp);
490   safr = std::fabs(dRho*salf);                    609   safr = std::fabs(dRho*salf);
491   if(safr>safz) { safz=safr; }                    610   if(safr>safz) { safz=safr; }
492                                                   611 
493   return safz;                                    612   return safz;
494 }                                                 613 }
495                                                   614 
496 //////////////////////////////////////////////    615 ///////////////////////////////////////////////////////////////////////////////
497 //                                                616 //
498 // Calculate distance to surface of shape from    617 // Calculate distance to surface of shape from 'inside'
499 //                                             << 618 
500 G4double G4Paraboloid::DistanceToOut(const G4T    619 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p,
501                                     const G4Th    620                                     const G4ThreeVector& v,
502                                     const G4bo    621                                     const G4bool calcNorm,
503                                           G4bo << 622                                           G4bool *validNorm,
504                                           G4Th << 623                                           G4ThreeVector *n  ) const
505 {                                                 624 {
506   G4double rho2 = p.perp2(), paraRho2 = std::f    625   G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
507   G4double vRho2 = v.perp2(), intersection;       626   G4double vRho2 = v.perp2(), intersection;
508   G4double tol2 = kCarTolerance*kCarTolerance;    627   G4double tol2 = kCarTolerance*kCarTolerance;
509   G4double tolh = 0.5*kCarTolerance;              628   G4double tolh = 0.5*kCarTolerance;
510                                                   629 
511   if(calcNorm) { *validNorm = false; }            630   if(calcNorm) { *validNorm = false; }
512                                                   631 
513   // We have that the particle p follows the l    632   // We have that the particle p follows the line x = p + s * v
514   // meaning x = p.x() + s * v.x(), y = p.y()     633   // meaning x = p.x() + s * v.x(), y = p.y() + s * v.y() and
515   // z = p.z() + s * v.z()                        634   // z = p.z() + s * v.z()
516   // The equation for all points on the surfac    635   // The equation for all points on the surface (surface expanded for
517   // to include all z) x^2 + y^2 = k1 * z + k2    636   // to include all z) x^2 + y^2 = k1 * z + k2 => .. =>
518   // => s = (A +- std::sqrt(A^2 + B)) / vRho2     637   // => s = (A +- std::sqrt(A^2 + B)) / vRho2
519   // where:                                       638   // where:
520   //                                              639   //
521   G4double A = k1 / 2 * v.z() - p.x() * v.x()     640   G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y();
522   //                                              641   //
523   // and:                                         642   // and:
524   //                                              643   //
525   G4double B = (-rho2 + paraRho2) * vRho2;        644   G4double B = (-rho2 + paraRho2) * vRho2;
526                                                   645 
527   if ( rho2 < paraRho2 && sqr(rho2 - paraRho2     646   if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2
528     && std::fabs(p.z()) < dz - kCarTolerance)     647     && std::fabs(p.z()) < dz - kCarTolerance)
529   {                                               648   {
530     // Make sure it's safely inside.              649     // Make sure it's safely inside.
531                                                   650 
532     if(v.z() > 0)                                 651     if(v.z() > 0)
533     {                                             652     {
534       // It's heading upwards, check where it     653       // It's heading upwards, check where it colides with the plane z = dz.
535       // When it does, is that in the surface     654       // When it does, is that in the surface of the paraboloid.
536       // z = p.z() + variable * v.z() for all     655       // z = p.z() + variable * v.z() for all points the particle can go.
537       // => variable = (z - p.z()) / v.z() so     656       // => variable = (z - p.z()) / v.z() so intersection must be:
538                                                   657 
539       intersection = (dz - p.z()) / v.z();        658       intersection = (dz - p.z()) / v.z();
540       G4ThreeVector ip = p + intersection * v;    659       G4ThreeVector ip = p + intersection * v; // Point of intersection.
541                                                   660 
542       if(ip.perp2() < sqr(r2 + kCarTolerance))    661       if(ip.perp2() < sqr(r2 + kCarTolerance))
543       {                                           662       {
544         if(calcNorm)                              663         if(calcNorm)
545         {                                         664         {
546           *n = G4ThreeVector(0, 0, 1);            665           *n = G4ThreeVector(0, 0, 1);
547           if(r2 < tolh || ip.perp2() > sqr(r2     666           if(r2 < tolh || ip.perp2() > sqr(r2 - tolh))
548           {                                       667           {
549             *n += G4ThreeVector(ip.x(), ip.y()    668             *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
550             *n = n->unit();                       669             *n = n->unit();
551           }                                       670           }
552           *validNorm = true;                      671           *validNorm = true;
553         }                                         672         }
554         return intersection;                      673         return intersection;
555       }                                           674       }
556     }                                             675     }
557     else if(v.z() < 0)                            676     else if(v.z() < 0)
558     {                                             677     {
559       // It's heading downwards, check were it    678       // It's heading downwards, check were it colides with the plane z = -dz.
560       // When it does, is that in the surface     679       // When it does, is that in the surface of the paraboloid.
561       // z = p.z() + variable * v.z() for all     680       // z = p.z() + variable * v.z() for all points the particle can go.
562       // => variable = (z - p.z()) / v.z() so     681       // => variable = (z - p.z()) / v.z() so intersection must be:
563                                                   682 
564       intersection = (-dz - p.z()) / v.z();       683       intersection = (-dz - p.z()) / v.z();
565       G4ThreeVector ip = p + intersection * v;    684       G4ThreeVector ip = p + intersection * v; // Point of intersection.
566                                                   685 
567       if(ip.perp2() < sqr(r1 + tolh))             686       if(ip.perp2() < sqr(r1 + tolh))
568       {                                           687       {
569         if(calcNorm)                              688         if(calcNorm)
570         {                                         689         {
571           *n = G4ThreeVector(0, 0, -1);           690           *n = G4ThreeVector(0, 0, -1);
572           if(r1 < tolh || ip.perp2() > sqr(r1     691           if(r1 < tolh || ip.perp2() > sqr(r1 - tolh))
573           {                                       692           {
574             *n += G4ThreeVector(ip.x(), ip.y()    693             *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
575             *n = n->unit();                       694             *n = n->unit();
576           }                                       695           }
577           *validNorm = true;                      696           *validNorm = true;
578         }                                         697         }
579         return intersection;                      698         return intersection;
580       }                                           699       }
581     }                                             700     }
582                                                   701 
583     // Now check for collisions with paraboloi    702     // Now check for collisions with paraboloid surface.
584                                                   703 
585     if(vRho2 == 0) // Needs to be treated sepe    704     if(vRho2 == 0) // Needs to be treated seperately.
586     {                                             705     {
587       intersection = ((rho2 - k2)/k1 - p.z())/    706       intersection = ((rho2 - k2)/k1 - p.z())/v.z();
588       if(calcNorm)                                707       if(calcNorm)
589       {                                           708       {
590         G4ThreeVector intersectionP = p + v *     709         G4ThreeVector intersectionP = p + v * intersection;
591         *n = G4ThreeVector(intersectionP.x(),     710         *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
592         *n = n->unit();                           711         *n = n->unit();
593                                                   712 
594         *validNorm = true;                        713         *validNorm = true;
595       }                                           714       }
596       return intersection;                        715       return intersection;
597     }                                             716     }
598     else if( ((A <= 0) && (B >= sqr(A) * (sqr(    717     else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0))
599     {                                             718     {
600       // intersection = (A + std::sqrt(B + sqr    719       // intersection = (A + std::sqrt(B + sqr(A))) / vRho2;
601       // The above calculation has a precision    720       // The above calculation has a precision problem:
602       // known problem of solving quadratic eq    721       // known problem of solving quadratic equation with small A  
603                                                   722 
604       A = A/vRho2;                                723       A = A/vRho2;
605       B = (k1 * p.z() + k2 - rho2)/vRho2;         724       B = (k1 * p.z() + k2 - rho2)/vRho2;
606       intersection = B/(-A + std::sqrt(B + sqr    725       intersection = B/(-A + std::sqrt(B + sqr(A)));
607       if(calcNorm)                                726       if(calcNorm)
608       {                                           727       {
609         G4ThreeVector intersectionP = p + v *     728         G4ThreeVector intersectionP = p + v * intersection;
610         *n = G4ThreeVector(intersectionP.x(),     729         *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
611         *n = n->unit();                           730         *n = n->unit();
612         *validNorm = true;                        731         *validNorm = true;
613       }                                           732       }
614       return intersection;                        733       return intersection;
615     }                                             734     }
616     std::ostringstream message;                   735     std::ostringstream message;
617     message << "There is no intersection betwe    736     message << "There is no intersection between given line and solid!"
618             << G4endl                             737             << G4endl
619             << "          p = " << p << G4endl    738             << "          p = " << p << G4endl
620             << "          v = " << v;             739             << "          v = " << v;
621     G4Exception("G4Paraboloid::DistanceToOut(p    740     G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
622                 JustWarning, message);            741                 JustWarning, message);
623                                                   742 
624     return kInfinity;                             743     return kInfinity;
625   }                                               744   }
626   else if ( (rho2 < paraRho2 + kCarTolerance      745   else if ( (rho2 < paraRho2 + kCarTolerance
627          || sqr(rho2 - paraRho2 - 0.25 * tol2)    746          || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 )
628          && std::fabs(p.z()) < dz + tolh)         747          && std::fabs(p.z()) < dz + tolh) 
629   {                                               748   {
630     // If this is true we're somewhere in the     749     // If this is true we're somewhere in the border.
631                                                   750     
632     G4ThreeVector normal = G4ThreeVector (p.x(    751     G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2);
633                                                   752 
634     if(std::fabs(p.z()) > dz - tolh)              753     if(std::fabs(p.z()) > dz - tolh)
635     {                                             754     {
636       // We're in the lower or upper edge         755       // We're in the lower or upper edge
637       //                                          756       //
638       if( ((v.z() > 0) && (p.z() > 0)) || ((v.    757       if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) )
639       {             // If we're heading out of    758       {             // If we're heading out of the object that is treated here
640         if(calcNorm)                              759         if(calcNorm)
641         {                                         760         {
642           *validNorm = true;                      761           *validNorm = true;
643           if(p.z() > 0)                           762           if(p.z() > 0)
644             { *n = G4ThreeVector(0, 0, 1); }      763             { *n = G4ThreeVector(0, 0, 1); }
645           else                                    764           else
646             { *n = G4ThreeVector(0, 0, -1); }     765             { *n = G4ThreeVector(0, 0, -1); }
647         }                                         766         }
648         return 0;                                 767         return 0;
649       }                                           768       }
650                                                   769 
651       if(v.z() == 0)                              770       if(v.z() == 0)
652       {                                           771       {
653         // Case where we're moving inside the     772         // Case where we're moving inside the surface needs to be
654         // treated separately.                    773         // treated separately.
655         // Distance until it goes out through     774         // Distance until it goes out through a side is returned.
656                                                   775 
657         G4double r = (p.z() > 0)? r2 : r1;        776         G4double r = (p.z() > 0)? r2 : r1;
658         G4double pDotV = p.dot(v);                777         G4double pDotV = p.dot(v);
659         A = vRho2 * ( sqr(r) - sqr(p.x()) - sq    778         A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y()));
660         intersection = (-pDotV + std::sqrt(A +    779         intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2;
661                                                   780 
662         if(calcNorm)                              781         if(calcNorm)
663         {                                         782         {
664           *validNorm = true;                      783           *validNorm = true;
665                                                   784 
666           *n = (G4ThreeVector(0, 0, p.z()/std:    785           *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z()))
667               + G4ThreeVector(p.x() + v.x() *     786               + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y()
668               * intersection, -k1/2).unit()).u    787               * intersection, -k1/2).unit()).unit();
669         }                                         788         }
670         return intersection;                      789         return intersection;
671       }                                           790       }
672     }                                             791     }
673     //                                            792     //
674     // Problem in the Logic :: Following condi    793     // Problem in the Logic :: Following condition for point on upper surface 
675     //                         and Vz<0  will     794     //                         and Vz<0  will return 0 (Problem #1015), but
676     //                         it has to retur    795     //                         it has to return intersection with parabolic
677     //                         surface or with    796     //                         surface or with lower plane surface (z = -dz)
678     // The logic has to be  :: If not found in    797     // The logic has to be  :: If not found intersection until now,
679     // do not exit but continue to search for     798     // do not exit but continue to search for possible intersection.
680     // Only for point situated on both borders    799     // Only for point situated on both borders (Z and parabolic)
681     // this condition has to be taken into acc    800     // this condition has to be taken into account and done later
682     //                                            801     //
683     //                                            802     //
684     // else if(normal.dot(v) >= 0)                803     // else if(normal.dot(v) >= 0)
685     // {                                          804     // {
686     //   if(calcNorm)                             805     //   if(calcNorm)
687     //   {                                        806     //   {
688     //     *validNorm = true;                     807     //     *validNorm = true;
689     //     *n = normal.unit();                    808     //     *n = normal.unit();
690     //   }                                        809     //   }
691     //   return 0;                                810     //   return 0;
692     // }                                          811     // }
693                                                   812 
694     if(v.z() > 0)                                 813     if(v.z() > 0)
695     {                                             814     {
696       // Check for collision with upper edge.     815       // Check for collision with upper edge.
697                                                   816 
698       intersection = (dz - p.z()) / v.z();        817       intersection = (dz - p.z()) / v.z();
699       G4ThreeVector ip = p + intersection * v;    818       G4ThreeVector ip = p + intersection * v;
700                                                   819 
701       if(ip.perp2() < sqr(r2 - tolh))             820       if(ip.perp2() < sqr(r2 - tolh))
702       {                                           821       {
703         if(calcNorm)                              822         if(calcNorm)
704         {                                         823         {
705           *validNorm = true;                      824           *validNorm = true;
706           *n = G4ThreeVector(0, 0, 1);            825           *n = G4ThreeVector(0, 0, 1);
707         }                                         826         }
708         return intersection;                      827         return intersection;
709       }                                           828       }
710       else if(ip.perp2() < sqr(r2 + tolh))        829       else if(ip.perp2() < sqr(r2 + tolh))
711       {                                           830       {
712         if(calcNorm)                              831         if(calcNorm)
713         {                                         832         {
714           *validNorm = true;                      833           *validNorm = true;
715           *n = G4ThreeVector(0, 0, 1)             834           *n = G4ThreeVector(0, 0, 1)
716              + G4ThreeVector(ip.x(), ip.y(), -    835              + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
717           *n = n->unit();                         836           *n = n->unit();
718         }                                         837         }
719         return intersection;                      838         return intersection;
720       }                                           839       }
721     }                                             840     }
722     if( v.z() < 0)                                841     if( v.z() < 0)
723     {                                             842     {
724       // Check for collision with lower edge.     843       // Check for collision with lower edge.
725                                                   844 
726       intersection = (-dz - p.z()) / v.z();       845       intersection = (-dz - p.z()) / v.z();
727       G4ThreeVector ip = p + intersection * v;    846       G4ThreeVector ip = p + intersection * v;
728                                                   847 
729       if(ip.perp2() < sqr(r1 - tolh))             848       if(ip.perp2() < sqr(r1 - tolh))
730       {                                           849       {
731         if(calcNorm)                              850         if(calcNorm)
732         {                                         851         {
733           *validNorm = true;                      852           *validNorm = true;
734           *n = G4ThreeVector(0, 0, -1);           853           *n = G4ThreeVector(0, 0, -1);
735         }                                         854         }
736         return intersection;                      855         return intersection;
737       }                                           856       }
738       else if(ip.perp2() < sqr(r1 + tolh))        857       else if(ip.perp2() < sqr(r1 + tolh))
739       {                                           858       {
740         if(calcNorm)                              859         if(calcNorm)
741         {                                         860         {
742           *validNorm = true;                      861           *validNorm = true;
743           *n = G4ThreeVector(0, 0, -1)            862           *n = G4ThreeVector(0, 0, -1)
744              + G4ThreeVector(ip.x(), ip.y(), -    863              + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
745           *n = n->unit();                         864           *n = n->unit();
746         }                                         865         }
747         return intersection;                      866         return intersection;
748       }                                           867       }
749     }                                             868     }
750                                                   869 
751     // Note: comparison with zero below would     870     // Note: comparison with zero below would not be correct !
752     //                                            871     //
753     if(std::fabs(vRho2) > tol2) // precision e    872     if(std::fabs(vRho2) > tol2) // precision error in the calculation of
754     {                           // intersectio    873     {                           // intersection = (A+std::sqrt(B+sqr(A)))/vRho2
755       A = A/vRho2;                                874       A = A/vRho2;
756       B = (k1 * p.z() + k2 - rho2);               875       B = (k1 * p.z() + k2 - rho2);
757       if(std::fabs(B)>kCarTolerance)              876       if(std::fabs(B)>kCarTolerance)
758       {                                           877       {
759         B = (B)/vRho2;                            878         B = (B)/vRho2;
760         intersection = B/(-A + std::sqrt(B + s    879         intersection = B/(-A + std::sqrt(B + sqr(A)));
761       }                                           880       }
762       else                      // Point is On    881       else                      // Point is On both borders: Z and parabolic
763       {                         // solution de    882       {                         // solution depends on normal.dot(v) sign
764         if(normal.dot(v) >= 0)                    883         if(normal.dot(v) >= 0)
765         {                                         884         {
766           if(calcNorm)                            885           if(calcNorm)
767           {                                       886           {
768             *validNorm = true;                    887             *validNorm = true;
769             *n = normal.unit();                   888             *n = normal.unit();
770           }                                       889           }
771           return 0;                               890           return 0;
772         }                                         891         }
773         intersection = 2.*A;                      892         intersection = 2.*A;
774       }                                           893       }
775     }                                             894     }
776     else                                          895     else
777     {                                             896     {
778       intersection = ((rho2 - k2) / k1 - p.z()    897       intersection = ((rho2 - k2) / k1 - p.z()) / v.z();
779     }                                             898     }
780                                                   899 
781     if(calcNorm)                                  900     if(calcNorm)
782     {                                             901     {
783       *validNorm = true;                          902       *validNorm = true;
784       *n = G4ThreeVector(p.x() + intersection     903       *n = G4ThreeVector(p.x() + intersection * v.x(), p.y()
785          + intersection * v.y(), - k1 / 2);       904          + intersection * v.y(), - k1 / 2);
786       *n = n->unit();                             905       *n = n->unit();
787     }                                             906     }
788     return intersection;                          907     return intersection;
789   }                                               908   }
790   else                                            909   else
791   {                                               910   {
792 #ifdef G4SPECSDEBUG                               911 #ifdef G4SPECSDEBUG
793     if(kOutside == Inside(p))                     912     if(kOutside == Inside(p))
794     {                                             913     {
795       G4Exception("G4Paraboloid::DistanceToOut    914       G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
796                   JustWarning, "Point p is out    915                   JustWarning, "Point p is outside!");
797     }                                             916     }
798     else                                          917     else
799       G4Exception("G4Paraboloid::DistanceToOut    918       G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
800                   JustWarning, "There's an err    919                   JustWarning, "There's an error in this functions code.");
801 #endif                                            920 #endif
802     return kInfinity;                             921     return kInfinity;
803   }                                               922   }
804   return 0;                                       923   return 0;
805 }                                                 924 } 
806                                                   925 
807 //////////////////////////////////////////////    926 ///////////////////////////////////////////////////////////////////////////////
808 //                                                927 //
809 // Calculate distance (<=actual) to closest su    928 // Calculate distance (<=actual) to closest surface of shape from inside
810 //                                             << 929 
811 G4double G4Paraboloid::DistanceToOut(const G4T    930 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p) const
812 {                                                 931 {
813   G4double safe=0.0,rho,safeR,safeZ ;             932   G4double safe=0.0,rho,safeR,safeZ ;
814   G4double tanRMax,secRMax,pRMax ;                933   G4double tanRMax,secRMax,pRMax ;
815                                                   934 
816 #ifdef G4SPECSDEBUG                               935 #ifdef G4SPECSDEBUG
817   if( Inside(p) == kOutside )                     936   if( Inside(p) == kOutside )
818   {                                               937   {
819      G4cout << G4endl ;                           938      G4cout << G4endl ;
820      DumpInfo();                                  939      DumpInfo();
821      std::ostringstream message;                  940      std::ostringstream message;
822      G4long oldprc = message.precision(16);    << 941      G4int oldprc = message.precision(16);
823      message << "Point p is outside !?" << G4e    942      message << "Point p is outside !?" << G4endl
824              << "Position:" << G4endl             943              << "Position:" << G4endl
825              << "   p.x() = "   << p.x()/mm <<    944              << "   p.x() = "   << p.x()/mm << " mm" << G4endl
826              << "   p.y() = "   << p.y()/mm <<    945              << "   p.y() = "   << p.y()/mm << " mm" << G4endl
827              << "   p.z() = "   << p.z()/mm <<    946              << "   p.z() = "   << p.z()/mm << " mm";
828      message.precision(oldprc) ;                  947      message.precision(oldprc) ;
829      G4Exception("G4Paraboloid::DistanceToOut(    948      G4Exception("G4Paraboloid::DistanceToOut(p)", "GeomSolids1002",
830                  JustWarning, message);           949                  JustWarning, message);
831   }                                               950   }
832 #endif                                            951 #endif
833                                                   952 
834   rho = p.perp();                                 953   rho = p.perp();
835   safeZ = dz - std::fabs(p.z()) ;                 954   safeZ = dz - std::fabs(p.z()) ;
836                                                   955 
837   tanRMax = (r2 - r1)*0.5/dz ;                    956   tanRMax = (r2 - r1)*0.5/dz ;
838   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;    957   secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
839   pRMax   = tanRMax*p.z() + (r1+r2)*0.5 ;         958   pRMax   = tanRMax*p.z() + (r1+r2)*0.5 ;
840   safeR  = (pRMax - rho)/secRMax ;                959   safeR  = (pRMax - rho)/secRMax ;
841                                                   960 
842   if (safeZ < safeR) { safe = safeZ; }            961   if (safeZ < safeR) { safe = safeZ; }
843   else { safe = safeR; }                          962   else { safe = safeR; }
844   if ( safe < 0.5 * kCarTolerance ) { safe = 0    963   if ( safe < 0.5 * kCarTolerance ) { safe = 0; }
845   return safe ;                                   964   return safe ;
846 }                                                 965 }
847                                                   966 
848 //////////////////////////////////////////////    967 //////////////////////////////////////////////////////////////////////////
849 //                                                968 //
850 // G4EntityType                                   969 // G4EntityType
851 //                                             << 970 
852 G4GeometryType G4Paraboloid::GetEntityType() c    971 G4GeometryType G4Paraboloid::GetEntityType() const
853 {                                                 972 {
854   return {"G4Paraboloid"};                     << 973   return G4String("G4Paraboloid");
855 }                                                 974 }
856                                                   975 
857 //////////////////////////////////////////////    976 //////////////////////////////////////////////////////////////////////////
858 //                                                977 //
859 // Make a clone of the object                     978 // Make a clone of the object
860 //                                             << 979 
861 G4VSolid* G4Paraboloid::Clone() const             980 G4VSolid* G4Paraboloid::Clone() const
862 {                                                 981 {
863   return new G4Paraboloid(*this);                 982   return new G4Paraboloid(*this);
864 }                                                 983 }
865                                                   984 
866 //////////////////////////////////////////////    985 //////////////////////////////////////////////////////////////////////////
867 //                                                986 //
868 // Stream object contents to an output stream     987 // Stream object contents to an output stream
869 //                                             << 988 
870 std::ostream& G4Paraboloid::StreamInfo( std::o    989 std::ostream& G4Paraboloid::StreamInfo( std::ostream& os ) const
871 {                                                 990 {
872   G4long oldprc = os.precision(16);            << 991   G4int oldprc = os.precision(16);
873   os << "-------------------------------------    992   os << "-----------------------------------------------------------\n"
874      << "    *** Dump for solid - " << GetName    993      << "    *** Dump for solid - " << GetName() << " ***\n"
875      << "    =================================    994      << "    ===================================================\n"
876      << " Solid type: G4Paraboloid\n"             995      << " Solid type: G4Paraboloid\n"
877      << " Parameters: \n"                         996      << " Parameters: \n"
878      << "    z half-axis:   " << dz/mm << " mm    997      << "    z half-axis:   " << dz/mm << " mm \n"
879      << "    radius at -dz: " << r1/mm << " mm    998      << "    radius at -dz: " << r1/mm << " mm \n"
880      << "    radius at dz:  " << r2/mm << " mm    999      << "    radius at dz:  " << r2/mm << " mm \n"
881      << "-------------------------------------    1000      << "-----------------------------------------------------------\n";
882   os.precision(oldprc);                           1001   os.precision(oldprc);
883                                                   1002 
884   return os;                                      1003   return os;
885 }                                                 1004 }
886                                                   1005 
887 //////////////////////////////////////////////    1006 ////////////////////////////////////////////////////////////////////
888 //                                                1007 //
889 // GetPointOnSurface                              1008 // GetPointOnSurface
890 //                                             << 1009 
891 G4ThreeVector G4Paraboloid::GetPointOnSurface(    1010 G4ThreeVector G4Paraboloid::GetPointOnSurface() const
892 {                                                 1011 {
893   G4double A = (fSurfaceArea == 0)? CalculateS    1012   G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea;
894   G4double z = G4RandFlat::shoot(0.,1.);       << 1013   G4double z = RandFlat::shoot(0.,1.);
895   G4double phi = G4RandFlat::shoot(0., twopi); << 1014   G4double phi = RandFlat::shoot(0., twopi);
896   if(pi*(sqr(r1) + sqr(r2))/A >= z)               1015   if(pi*(sqr(r1) + sqr(r2))/A >= z)
897   {                                               1016   {
898     G4double rho;                                 1017     G4double rho;
899     if(pi * sqr(r1) / A > z)                      1018     if(pi * sqr(r1) / A > z)
900     {                                             1019     {
901       rho = r1 * std::sqrt(G4RandFlat::shoot(0 << 1020       rho = r1 * std::sqrt(RandFlat::shoot(0., 1.));
902       return { rho * std::cos(phi), rho * std: << 1021       return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz);
903     }                                             1022     }
904     else                                          1023     else
905     {                                             1024     {
906       rho = r2 * std::sqrt(G4RandFlat::shoot(0 << 1025       rho = r2 * std::sqrt(RandFlat::shoot(0., 1));
907       return { rho * std::cos(phi), rho * std: << 1026       return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), dz);
908     }                                             1027     }
909   }                                               1028   }
910   else                                            1029   else
911   {                                               1030   {
912     z = G4RandFlat::shoot(0., 1.)*2*dz - dz;   << 1031     z = RandFlat::shoot(0., 1.)*2*dz - dz;
913     return { std::sqrt(z*k1 + k2)*std::cos(phi << 1032     return G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi),
914              std::sqrt(z*k1 + k2)*std::sin(phi << 1033                          std::sqrt(z*k1 + k2)*std::sin(phi), z);
                                                   >> 1034   }
                                                   >> 1035 }
                                                   >> 1036 
                                                   >> 1037 G4ThreeVectorList*
                                                   >> 1038 G4Paraboloid::CreateRotatedVertices(const G4AffineTransform& pTransform,
                                                   >> 1039                                           G4int& noPolygonVertices) const
                                                   >> 1040 {
                                                   >> 1041   G4ThreeVectorList *vertices;
                                                   >> 1042   G4ThreeVector vertex;
                                                   >> 1043   G4double meshAnglePhi, cosMeshAnglePhiPer2,
                                                   >> 1044            crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi,
                                                   >> 1045            sRho, dRho, rho, lastRho = 0., swapRho;
                                                   >> 1046   G4double rx, ry, rz, k3, k4, zm;
                                                   >> 1047   G4int crossSectionPhi, noPhiCrossSections, noRhoSections;
                                                   >> 1048 
                                                   >> 1049   // Phi cross sections
                                                   >> 1050   //
                                                   >> 1051   noPhiCrossSections = G4int(twopi/kMeshAngleDefault)+1;  // =9!
                                                   >> 1052 /*
                                                   >> 1053   if (noPhiCrossSections<kMinMeshSections)          // <3
                                                   >> 1054   {
                                                   >> 1055     noPhiCrossSections=kMinMeshSections;
                                                   >> 1056   }
                                                   >> 1057   else if (noPhiCrossSections>kMaxMeshSections)     // >37
                                                   >> 1058   {
                                                   >> 1059     noPhiCrossSections=kMaxMeshSections;
                                                   >> 1060   }
                                                   >> 1061 */
                                                   >> 1062   meshAnglePhi=twopi/(noPhiCrossSections-1);
                                                   >> 1063 
                                                   >> 1064   sAnglePhi = -meshAnglePhi*0.5*0;
                                                   >> 1065   cosMeshAnglePhiPer2 = std::cos(meshAnglePhi / 2.);
                                                   >> 1066 
                                                   >> 1067   noRhoSections = G4int(pi/2/kMeshAngleDefault) + 1;
                                                   >> 1068 
                                                   >> 1069   // There is no obvious value for noRhoSections, at the moment the parabola is
                                                   >> 1070   // viewed as a quarter circle mean this formula for it.
                                                   >> 1071 
                                                   >> 1072   // An alternetive would be to calculate max deviation from parabola and
                                                   >> 1073   // keep adding new vertices there until it was under a decided constant.
                                                   >> 1074 
                                                   >> 1075   // maxDeviation on a line between points (rho1, z1) and (rho2, z2) is given
                                                   >> 1076   // by rhoMax = sqrt(k1 * z + k2) - z * (rho2 - rho1)
                                                   >> 1077   //           / (z2 - z1) - (rho1 * z2 - rho2 * z1) / (z2 - z1)
                                                   >> 1078   // where z is k1 / 2 * (rho1 + rho2) - k2 / k1
                                                   >> 1079 
                                                   >> 1080   sRho = r1;
                                                   >> 1081   dRho = (r2 - r1) / double(noRhoSections - 1);
                                                   >> 1082 
                                                   >> 1083   vertices=new G4ThreeVectorList();
                                                   >> 1084 
                                                   >> 1085   if (vertices)
                                                   >> 1086   {
                                                   >> 1087     for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
                                                   >> 1088          crossSectionPhi++)
                                                   >> 1089     {
                                                   >> 1090       crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
                                                   >> 1091       coscrossAnglePhi=std::cos(crossAnglePhi);
                                                   >> 1092       sincrossAnglePhi=std::sin(crossAnglePhi);
                                                   >> 1093       lastRho = 0;
                                                   >> 1094       for (int iRho=0; iRho < noRhoSections;
                                                   >> 1095            iRho++)
                                                   >> 1096       {
                                                   >> 1097         // Compute coordinates of cross section at section crossSectionPhi
                                                   >> 1098         //
                                                   >> 1099         if(iRho == noRhoSections - 1)
                                                   >> 1100         {
                                                   >> 1101           rho = r2;
                                                   >> 1102         }
                                                   >> 1103         else
                                                   >> 1104         {
                                                   >> 1105           rho = iRho * dRho + sRho;
                                                   >> 1106 
                                                   >> 1107           // This part is to ensure that the vertices
                                                   >> 1108           // will form a volume larger than the paraboloid
                                                   >> 1109 
                                                   >> 1110           k3 = k1 / (2*rho + dRho);
                                                   >> 1111           k4 = rho - k3 * (sqr(rho) - k2) / k1;
                                                   >> 1112           zm = (sqr(k1 / (2 * k3)) - k2) / k1;
                                                   >> 1113           rho += std::sqrt(k1 * zm + k2) - zm * k3 - k4;
                                                   >> 1114         }
                                                   >> 1115 
                                                   >> 1116         rho += (1 / cosMeshAnglePhiPer2 - 1) * (iRho * dRho + sRho);
                                                   >> 1117 
                                                   >> 1118         if(rho < lastRho)
                                                   >> 1119         {
                                                   >> 1120           swapRho = lastRho;
                                                   >> 1121           lastRho = rho + dRho;
                                                   >> 1122           rho = swapRho;
                                                   >> 1123         }
                                                   >> 1124         else
                                                   >> 1125         {
                                                   >> 1126           lastRho = rho + dRho;
                                                   >> 1127         }
                                                   >> 1128 
                                                   >> 1129         rx = coscrossAnglePhi*rho;
                                                   >> 1130         ry = sincrossAnglePhi*rho;
                                                   >> 1131         rz = (sqr(iRho * dRho + sRho) - k2) / k1;
                                                   >> 1132         vertex = G4ThreeVector(rx,ry,rz);
                                                   >> 1133         vertices->push_back(pTransform.TransformPoint(vertex));
                                                   >> 1134       }
                                                   >> 1135     }    // Phi
                                                   >> 1136     noPolygonVertices = noRhoSections ;
                                                   >> 1137   }
                                                   >> 1138   else
                                                   >> 1139   {
                                                   >> 1140     DumpInfo();
                                                   >> 1141     G4Exception("G4Paraboloid::CreateRotatedVertices()",
                                                   >> 1142                 "GeomSolids0003", FatalException,
                                                   >> 1143                 "Error in allocation of vertices. Out of memory !");
915   }                                               1144   }
                                                   >> 1145   return vertices;
916 }                                                 1146 }
917                                                   1147 
918 //////////////////////////////////////////////    1148 /////////////////////////////////////////////////////////////////////////////
919 //                                                1149 //
920 // Methods for visualisation                      1150 // Methods for visualisation
921 //                                             << 1151 
922 void G4Paraboloid::DescribeYourselfTo (G4VGrap    1152 void G4Paraboloid::DescribeYourselfTo (G4VGraphicsScene& scene) const
923 {                                                 1153 {
924   scene.AddSolid(*this);                          1154   scene.AddSolid(*this);
925 }                                                 1155 }
926                                                   1156 
927 G4Polyhedron* G4Paraboloid::CreatePolyhedron (    1157 G4Polyhedron* G4Paraboloid::CreatePolyhedron () const
928 {                                                 1158 {
929   return new G4PolyhedronParaboloid(r1, r2, dz    1159   return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi);
930 }                                                 1160 }
931                                                   1161 
                                                   >> 1162 
932 G4Polyhedron* G4Paraboloid::GetPolyhedron () c    1163 G4Polyhedron* G4Paraboloid::GetPolyhedron () const
933 {                                                 1164 {
934   if (fpPolyhedron == nullptr ||               << 1165   if (!fpPolyhedron ||
935       fRebuildPolyhedron ||                       1166       fRebuildPolyhedron ||
936       fpPolyhedron->GetNumberOfRotationStepsAt    1167       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
937       fpPolyhedron->GetNumberOfRotationSteps()    1168       fpPolyhedron->GetNumberOfRotationSteps())
938   {                                               1169   {
939     G4AutoLock l(&polyhedronMutex);               1170     G4AutoLock l(&polyhedronMutex);
940     delete fpPolyhedron;                          1171     delete fpPolyhedron;
941     fpPolyhedron = CreatePolyhedron();            1172     fpPolyhedron = CreatePolyhedron();
942     fRebuildPolyhedron = false;                   1173     fRebuildPolyhedron = false;
943     l.unlock();                                   1174     l.unlock();
944   }                                               1175   }
945   return fpPolyhedron;                            1176   return fpPolyhedron;
946 }                                                 1177 }
947                                                << 
948 #endif // !defined(G4GEOM_USE_UPARABOLOID) ||  << 
949                                                   1178